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ABSTRACT 


A new theory and computer program for combustion instability analysis are 
presented herein. The basic theoretical foundation resides in the concept of 

entropy-controlled energy growth or decay. Third order perturbation expansion is 
performed on the entropy— controlled acoustic energy equation to obtain the first order 
integrodifferential equation for the energy growth factor in terms of the linear, second, and 
third order energy growth parameters. These parameters are calculated from 

Navier-Stokes solutions with time averages performed on as many Navier-Stokes time 
steps as required to cover at least one peak wave period. 

Applications are made for one-dimensional Navier-Stokes solution for the SSME 
thrust chamber with cross section area variations taken into account. It is shown that 
instability occurs when the mean pressure is raised to 2000 psi with 30% disturbances. 

Instability also arises when the mean pressure is set at 2935 psi with 20% disturbances. 

The system with mean pressures and disturbances more adverse than these cases has been 
shown to be unstable. 

The present theory has a great potential and all avenues of further studies will 


prove to be fruitful. 
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NOMENCLATURE 


c 


p 



1J 


Specific heat at constant pressure 

Body force vector 
Mass diffusivity 
Internal energy density 
Stagnation energy 
Body force 

Convective flux vector 

Dissipative vector 

Total enthalpy 

Pressure 
Gas constant 
Entropy 

Time dependent variable vector 
Velocity 

Mass fraction 

Energy growth rate parameters of first order, 

second order, and third order, respectively. 

Specific heat ratio 
Energy growth factor 
Thermal conductivity 
Viscosity 
Density 

Total stress tensor 
Viscous stress tensor 
Reaction rate 


Subscripts and Superscripts 
’ Fluctuation 

Time averaged mean quantity 
0 Reference state 
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1. INTRODUCTION 

Unstable waves may exhibit a linear behavior initially under the low mean pressure, but 
tend to oscillate nonlinearly as the mean pressure increases, resulting possibly in sawtooth 
wave forms. Multidimensional effects become significant as transverse modes contribute to 
instability. Chemical reactions, atomization, vaporization, and turbulent flow 
environments must also be considered. With these complications affecting the overall 
stability behavior, we come to the question: What is the most rigorous method of 
determining combustion instability? 

If time— dependent Navier-Stokes solutions for combustion capable of generating 
both linear and nonlinear wave oscillations are available, this information alone may 
provide qualitative interpretation of instability as to the tendency of possible energy 
growth or decay. However, they do not provide quantitative data for instability. Will 
there be, then, a "measure" of instability? In fact, there have been many attempts in 
seeking such data, the so-called "growth rate parameter" [1—5]. Unfortunately, they are 
normally limited to linear instability. 

In order to accommodate nonlinear behavior, multidimensionality, and complex 
flowfield phenomena, we introduce a new approach, the Entropy— Controlled— Instability 
(ECI) method. The concept is similar to Flandro [6] in which the energy balance method 
was used in deriving the expression for energy growth from the acoustic energy equation. 
The focal point of the present study is the entropy— controlled energy equation which 
automatically takes into account shock wave oscillations in determining energy growth for 
instability. The asymptotic perturbation expansions of all acoustic energy terms lead to 
the entropy— controlled— energy equation. Applying the Green— Gauss theorem and taking 
time averages, we derive the stability integrodifferential equation for the energy growth 
factor. This factor is solved in terms of growth rate parameters which are determined from 
the Navier-Stokes solution. 
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The advantage of the present method is to provide stability information during any time 
period of Navier-Stokes solutions. Stability prediction capability is, therefore, limited only 
by the Navier-Stokes solver. 

In the following, we shall describe the governing equations, derivation of stability 
integrodifferential equation, solution procedure, and one-dimensional example problems for 
validation of the theory. Extension to multi dimensions and more complex flow fields is 
achieved simply by adopting an appropriate Navier-Stokes solver. The present 
formulation of stability analysis remains unchanged. 


2. GOVERNING EQUATIONS 

2.1 Navier-Stokes Equations 

The most general conservation form of Navier Stokes equations is given by 
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and E is the stagnation energy 

E = e + jV iVi = c p T-£ + 
and f ki is the body force and qj is the heat flux vector. 

1j = - AT,j + p D k S H k Y k „ 

Here, A and D are the thermal conductivity and mass diffusivity, respectively. H k is the 
total enthalpy of species k, Y k is the mass fraction for the species k, and w k is the reaction 
rate for the species k. Example problems in this report do not include reacting flows. 

Solution of the Navier— Stokes equations is obtained using the Taylor— Galer kin 
finite element method. Details of the solution procedure are found in [7]. 


2.2 Entropy— Controlled Stability Equation 


Suppose that the Navier— Stokes solution has been obtained with the results exhibiting 
sawthooth waves. Our objective is to determine whether such waves are stable or unstable. 
To this end we examine the conservation form of the energy equation, 

|t (/> E ) + (/> Ev i “ *yVj),i = 0 (2) 

where the comma implies partial derivatives and is the stress tensor, 

*ij = - P^ij + M (vi.j + Vj,i - 1 v k , k ^j) (3) 

From thermodynamic relations it can be shown (appendix A) that 

P E >i = pP>i + f S >i + P v j v j>i ( 4 ) 

where S is the specific enthropy per unit mass. Substituting (4) into (2) yields 


|f 0> E ) + E (P v i),i + Vi + £ S,i + pv jVj)i — ( ^"ij v j ) , i = 0 


.(5) 


This is the entropy— controlled— energy equation, instrumental in determining the nonlinear 
instability. 
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Assuming that the Navier-Stokes solutions for density p, pressure p, and velocity Vj 
represent the sum of mean and fluctuation parts, we write 

p = p + p ( 6 ) 

P = P + P' CO 

vj = Vj + v- (8) 

where the symbols, bar and prime, denote the mean and perturbation quantities, 
respectively. 

From thermodynamic relations we may write the entropy difference in the form 


S - S 0 = Rin 




1 

I V - ^ -1 


or 


S = R 


S(D + S( 2) + S( 3) +S( 


4 ) 


+ s 0 (9) 

where S 0 represents the entropy at the initial state and S ( are given in Appendix B. 

Our objective is to establish quantitative criteria whether the system is stable or 
unstable when we are provided with the Navier-Stokes solution exhibiting wave 
oscillations during unsteady motions. To this end, let e be the energy growth factor, e > 0 
with e = 1 indicating the neutral stability. We then substitute (6) through (9) into (5), 
expand each term of the energy equation in terms of e, integrate by parts (or using 
Green— Gauss theorem), and take time averages 

Writing (5) in an integral form 

< J n [ft 0>E) + EOnr,),, + v, (E + | S„ + - (^Tj),,] dn = 0 


(10) 


Integrating (10) by parts, 

< J Q If (^ E ) dn + } r v i n i + v i (p n i + ft S n i + P v j v j n i) " a ij v j n i 


dr 


E ,i P V i + ( V i ?)»i + ( V i f).i S + (P V i V j)>i v j 


dn ) = o 


(ii) 
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where (•) implies time averages. A typical term in (11) for multiples of two or more 
variables appears in the form 



(•) dfi ) = ( (<$ 0 + 6(5, + e 2 6 2 + € 3 S 3 + ...) dO) 


(12) 


Here S 0 term contains only the mean quantity, S v the first order perturbation, <5 2 , the 
second order perturbation, etc. See detailed derivations in Appendix C. 

It follows from (12) that the perturbed acoustic equation takes the form 


Jf (t ! E, + f 3 E 2 + <%) = e ! I, + <>1, + e 4 I, (13) 

Thus, finally, the entropy— controlled stability equation becomes (See Appendix C) 

jjf- a 2 e2 “ a 3 f3 = 0 ( 14 ) 

where a v a 2 , and a 3 are growth rate parameters of first, second and third order, 
respectively. 
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I a = < f b<2)dQ ) — < [ ct^njdT) 
Jo Jr 

(17b) 

I, = ( [ b<3>dH ) - < [ c^JnjdT) 
Jo Jr 

(17c) 


where < • ) implies the time average and explicit forms of integrands are shown in Appendix 
D. It should be noted that all terms with Q represent acoustic energy in the domain 
whereas those with T denote acoustic intensities along the boundary surfaces. The linear 
growth rate parameter a x does not contain the terms associated with entropy whereas the 
nonlinear growth rate parameters a 2 and a 3 involve entropy— induced terms which are 
expected to play a role in energy dissipation leading to limit cycles and triggered 
instability. 

The basic ingredients of integrands in Eq. (15) are the data from Navier-Stokes 
solutions. The mean quantities are obtained as time averages of Navier-Stokes solutions 
within suitable time segments and the fluctuation (perturbation) quantities are the 
differences between the Navier-Stokes solutions and their time averages. 

To gain an insight into a solution of Eq. (14), we may neglect the last two terms of 
the left hand side of Eq. (14) and write 

(18) 

which yields a solution in the form 

bi £ = aqt + c t (19) 

To establish an initial condition we assume a neutral stability e = 1 at t = 0. This gives 
c t = 0. Thus, the solution takes the form 

e = eV . (20) 

* Under this initial condition, there exists a unique solution for any given aq with t > 0. It 
then follows that for stability we have 0 < e < 1 for — oo < cij < 0; for instability 1 < e < oo 
for 0 < a 1 < oo. 

i 
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Although these criteria are not applicable for the nonlinear equation (Eq. 14), similar 
initial conditions as postulated above can be used. That is, there exists a unique solution e 
for any given a t , a 2 , and a 3 with t > 0. 

Solutions of the nonlinear equation (Eq. 14) may be obtained using 
Newton-Raphson iterations. To this end, the residual of Eq. (14) is written as 


R- n + 


l>r 


_ At 

e n + l>r e n>r 2 


a l ( e n + l>r "b ^n>r) "b (^n+l’r "b e n>r) 


"b ^3 ( f n + l>r "b f n’r) 

The Newton-Raphson process for Eq. (24) takes the form 


( 21 ) 


^n + l’r ^ e n+l>r*l ^n+l’r 

where the Jacobian J n+1 , r becomes 

^ ^ | J j* ^ i 

+ = ^ — "5~" ( a t + ^**2 e n + L’r ^3 f n + l’r) 

C,t n + Pr 46 

and 

^^n + l’r + l — f n + l>r + 1 ~~ ^n + l>r 

Thus for each iterative step, we have 

^n + l’r + l ^n + l»r ^^n+l>r + l 


( 22 ) 


(23) 


(24) 


(25) 


The initial value for e begins with € n , r = 0 and £ n+l , r = 1. Iterations continue until 
convergence. 


3. SOLUTION PROCEDURE 

To solve the nonlinear ordinary differential equation (14), we proceed as follows: 


(1) With appropriate boundary and initial conditions, solve the Navier— Stokes 
equations using a numerical scheme capable of handling shock 
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discontinuities. Obtain p, Vj, and p. The Taylor-Galerkin Finite Element 
method is used in this study. 

(2) Advance time steps (At) of Navier-Stokes solutions to obtain wave 

oscillations to cover at least one wave period. 

(3) Take time averages for the period nAt (the range of n is approximately, 15 < 

n < 150, depending on frequencies f, n is small if f is high), corresponding to 

p, v it and p. 

(4) Calculate the fluctuation quantities as p' = p — p, v- = v i — v i; etc., where p, 
and Vj represent Navier-Stokes solutions. 

(5) Calculate the growth rate parameters a t , a 2 , and a 3 from (13a,b,c). 

(6) Solve the nonlinear ordinary differential equation (14) using the 

Newton— Rap hson method with a suitable initial guess for e. Ideally begin 
with e = 1, neutral stability. 

(7) Repeat steps 1 through 4 until the desired length of time has been advanced. 

Note that for each time— average period in step 4, above, instability and stability are 
determined by e > 1 and e < 1, respectively, with e = 1 being the neutral stability. If the 
system is found to be unstable, then it is not necessary to proceed to the next time step. 
However, for the entire ranges of time for which Navier-Stokes solutions are available, the 
stability analysis may be performed if desired, even if instability has been found in previous 
time steps. This is so because Navier-Stokes solutions are independent of the stability 
analysis as formulated here. Rather, the stability analysis in this formulation determines 
the state of stability or instability based on the current flowfield as calculated from the 
Navier-Stokes solution. 
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4. APPLICATIONS 

Our objective here is to prove validity of the present theory for combustion instability 
analysis. To this end, one dimensional nonreacting flow has been chosen for the geometry 
of SSME thrust chamber with cross section area variations taken into account (Fig. 1). 

The initial and boundary conditions for the Navier— Stokes solution consists of: 
Pressure p = p + d p sin (art + 0 o ) 

% disturbance, d = 10, 20, 30% 

mean pressure, p = 500, 2,000, 2,935 psi 

frequency, u) = 27rf > a/ 2L (L = distance between inlet and 

nozzle throat) 

Velocity (inlet) u = Ma (M = 0.2) 

Temperature T = 1000° R for p = 500 psi 

T = 4000° R for p = 2000 psi 
T = 6550° R for p = 2935 psi 
Other constants used in this analysis are: 

Specific heat ratio 7 = 1.2 

Mesh size Ax = 1.685 x 10‘ 2 m 

Courant Number C.N. = 0.6 

The computational time increment At is calculated at each time step of Navier-Stokes 
solution from the Courant number. Time averages for calculation of energy growth rate 
parameters a v a 2 , and a 3 are calculated over 15 to 120 intervals of Navier-Stokes At’s to 
cover at least an average of one peak at any grid point. For simplicity, viscosity is ignored 
in this example problem. The computer program listing is given in Appendix E. 

The Navier-Stokes solutions were obtained using the Taylor— Galerkin finite 
elements. Formulations of this method have been well documented and accuracy verified 
in the literature [7], 
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In Figs. 2 through 15, for each mean pressure and each % disturbance, the pressure and 
velocity oscillations are shown at various locations, x = — 0.31, 0.05 m, and 3.05 m, along 
with the corresponding energy growth factors versus time. 

In Figs. 2 and 3, (p = 500 psi, d = 10%, T = 1000° R), we note that it takes 0.0163 
sec. for the pressure at x = 3.06 m to begin decreasing and for the velocity to increase from 
zero. Notice that shock waves develop around t = 0.14 sec., but stability is maintained 
(e < 1) throughout since the mean pressure and disturbance are small. 

The response due to p = 500 psi, d = 30%, T = 1000° R, Figs. 4 and 5, is very 
similar to the case for d = 10%. Although the shock waves grow in magnitude and the 
energy growth factors increase, the system is still stable (e < 1). 

In Figs. 6 and 7, (p = 1000 psi, d = 20%, T = 4000° R), shock waves grow and the 
energy growth factors reach almost the level of neutral stability. But, instability has not 
been observed. 

The first instability has arrived at p = 2000 psi, d = 30%, T = 4000° R, Figs. 8 and 
9, in the time interval, 0.045 < t < 0.6 sec, where sawtooth type shock waves at x = 
3.06 m are prominent. 

In Figs. 10 and 11, in which the pressure is raised to p = 2935 psi with T = 6550° R, 
but disturbances are lowered to d = 10%, the system recovers stability. 

With the disturbances raised to d = 20%, p = 2935 psi, however, Figs. 12 and 13, 
notice that the energy growth factor rises sharply at t = 0.05 sec. where pressure decreases 
to almost zero, but shock waves rise rapidly. However, the instability for this case is not as 
severe as when pressure was lower (2000 psi) but disturbance was large (30%), as seen in 
^ Figs. 8 and 9. 

The most severe instability occurs when the disturbances are raised to d = 30%, 
with p = 2935 psi, Figs. 14 and 15. Notice that instability is spread over the wide time 
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range 0.045 < t < 0.06 sec., rather than a single peak for the case of d = 20% above. 
Similar situation existed for d = 30% with p = 2000 psi. It appears that instability is more 
sensitive to the increase in % disturbances than mean pressure. 

In Figs. 16 through 19, the energy growth rate parameters a v a 2 , and a 3 versus 
time are shown. When stable, the sum of a v a 2 , and a 3 is negative for the case of 
p = 500 psi and d = 10% (Fig. 16). If unstable, however, the sum is positive for cases of 
p = 2000 psi and d = 30% (Fig. 17), p = 2935 psi and d = 20% (Fig. 18), and p = 2935 psi 
and d = 30% (Fig. 19). Notice that as pressure increases the distribution of energy growth 
parameters become oscillatory. It is important to realize, however, that represents a 
linear instability whereas a 2 and a 3 contribute to the nonlinear instability as controlled by 
entropy. 

5. CONCLUSIONS 

To our knowledge, the full scale Navier— Stokes solutions combined with rigorous 
determinations of stability or instability during any time step of unsteady Navier-Stokes 
solutions have been carried out for the first time. The key to this success lies in the fact 
that the entropy is induced in the acoustic energy equation. It is shown that entropy is 
calculated automatically, contributing to the shock waves and instability. For small 
disturbances and low pressures the effect of entropy is negligible whereas it is activated 
freely when the mean pressures and disturbances are increased. 

To demonstrate the validity of the theory, the space shuttle main engine thrust 
chamber geometry was adopted for one dimensional flow but with cross section area 
variations taken into account. The computational results indicate that instability (e >.l) 
arises first when the mean pressure is raised to 2000 psi with 30% disturbances. Instability 
also arises when the mean pressure is set at 2935 psi with 20% disturbances. The system 
with mean pressures and disturbances more adverse than these quantities are shown to be 
unstable. 
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6. RECOMMENDATIONS 

Based on the studies reported herein the following recommendations are provided: 

(1) Extend the calculations to two-dimensional, axisymmetric cylindrical, and 
three dimensional geometries. 

(2) Investigate effects of chemical kinetics. 

(3) Investigate effects of Reynolds number (viscosity). 

(4) Investigate effects of atomization, vaporization, and spray droplet 
combustion. 

(5) Investigate effects of radiative heat transfer. 

In summary, it is the opinion of this principal investigator that the present theory 
has a great potential and all avenues of further studies will prove to be fruitful. 
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Fig. 5 Navier-Stokes solutions for velocity (a), (b) , (c) 
at various locations and (d) energy growth factors 
(£) versus time for p = 500 psi, d - 30%, T - 1000°R. 
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Fig. 6 Navier-Stokes solutions for pressure (a), (b), (c) 
at various locations and (d) energy growth factors 
(e) versus time for p * 2000 psi, d - 20Z, T * 4000°R 
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Fig. 7 Navier-Stokes solutions for velocity (a), (b) , (c) 
at various locations and (d) energy growth factors 
(s) versus time for p » 2000 psi, d - 20%, T * 4000°R 
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Fig. 9 Navier-Stokes solutions for velocity (a), (b) , (c) 
at various locations and (d) energy growth factors 
(e) versus time for p * 2000 psi, d * 30Z, T * 4000°R 
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Fig. 10 Navier-Stokes solutions for pressure (a), (b), (c) 
at various locations and (d) energy growth factors 
(e) versus time for p - 2935 psi, d - 10%, T - 6550°R. 
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Fig. 11 Navier-Stokes solutions for velocity (a), (b) , (c) 
at various locations and (d) energy growth factors 
(£) versus time for p * 2935 psi, d ■ 10Z, T * 6550°R. 
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Fig. 12 Navier-Stokes solutions for pressure (a), (b), (c) 
at various locations and (d) energy growth factors 
(e) versus time for p - 2935 psi, d - 20Z, T - 6550°R. 
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Fig. 3 


Navier-Stokes solutions for velocity (a), (b) , (c) 

at various locations and (d) energy growth factors 
(e) versus time for p = 500 psi, d ■ 10%, T = 1000°R. 
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Fig. 5 Navier-Stokes solutions for velocity (a), (b) , (c) 
at various locations and (d) energy growth factors 
(£) versus time for p = 500 psi, d = 30%, T = 1000°R. 
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7 Navier-Stokes solutions for velocity (a) , (b), (c) 
at various locations and (d) energy growth factors 
(s) versus time for p - 2000 psi, d * 20%, T * 4000°R 
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Fig. 9 Navier-Stokes solutions for velocity (a), (b) , (c) 
at various locations and (d) energy growth factors 
(e) versus time for p ■ 2000 psi, d ■ 30Z, T * 4000°R 
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Fig. 11 Navier-Stokes solutions for velocity (a), (b) , (c) 
at various locations and (d) energy growth factors 
(e) versus time for p ■ 2935 psi, d * 10%, T * 6550°R. 














Fig. 13 Navier-Stokes soluti 
at various locations 
(e) versus time for 

















TIME ( sec) 


Fig. 15 Navier-Stokes solutions for velocity (a), (b) , (c) 
at various locations and (d) energy growth factors 
(e) versus time for p - 2935 psi, d - 30%, T ■ 6550°R. 
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Fig. 16 Energy growth rate parameters a^, , ^3 

versus time, p = 500 psi, d = 10%, stable 
system (s<l), sum of a^, a^r a 3 less than 
zero . 
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Fig. 19 Energy growth rate parameters 
versus time, p = 2935 psi, d = 
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APPENDIX A 


DERIVATION OF ENERGY GRADIENTS IN TERMS OF ENTROPY GRADIENTS 

From an ideal gas law 

E- = (H.V exp & 

Po V</ K c v ' 


or 


(S-) = * (§- T + 


S-S r 


To 

Differentiating 


1 1 ( "K | 1 n 

n P>i = — {P )>i + — S >i 
P p ' L v 


or 


P Co 


P-i = c oP>i + ~T S ’i 
L P 

Now the gradient of the stagnation energy becomes 

E -i = ( c P T -f + J v j v j)>i 


or 


c v c v p 

E ’i + Ep P ’i - TTp /J ’i + V j V j Ji 
Substituting (A.l) into (A. 2), we obtain 

pE,i = p p , j + S,j + p v j v jii 


(A.l) 


(A.2) 


(A.3) 
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APPENDIX B 

Derivation of Entropy Perturbation 


S - S 0 = Rta 


(i + " (i + L) 

p t> 


\'K*1 


= R 


= R 


«ZT ^ (l + ?■) “ 61 (l + f) 

' p p 

' l n n n D ^ 


1 

p 


7-1 


5 

p p 


= *{(£*-£*) 

1 ' p ' p 


+ 5 


1 

55 


' p P J 


Thus 


S = R 


where 


S( 1) — 


S< 2) — 


s ( 3) - 5 


s c 1) + S( 2 ) + S (3 ) + S( 4) j + s 0 
1 El 1 £' 

L^p r*Jl 

&(*)*- £(f)*| 


L 7- 


£(*)*- 


S< 4 . - -Jg 


L 7- 



(B.l) 


(B.2) 


(B.3) 
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APPENDIX C 

DERIVATION OF INTEGRODIFFERENTIAL 
EQUATION FOR ENTROPY INDUCED ENERGY GROWTH 


From Eq. (11) and Eq. (A. 3) the energy equation takes the form 

|t (P E ) = " E (P v i)>i ~ v i [ p Pn + ft S >i + v j v j’i] + ( a ij v j)>i 

P v i P v i 

= — E(pv i ), i S fi — P, i - P v i v j v j>i + ( j V j ) > i 

1 1 V i 

= — (EpvJ,; + pvjEjj (pvjS),! + ^ S (pvjjj — — P, i 

PV jVj Vj , i + (^ijVj)^ 

= -(EpV^i-^pVjS)^ + (ffyVjJ.j 


+ 


1 P v i 

P v i E u + H S (P v i)>i ^-Pu-P v i v j v j>i 


Intetrating the above over the domain 0 and boundary T and taking the time averages 

P Vi 


r r r i pv^ 

< If(pE)dn) = <J [pv i E, i + K S(pv i ), i — P,{ ~ P v i v j v j>iJ dQ> 

J o J n L 

+ < - P v i E_ RP v i S + CTjjVj Sjdr) 

J r L 

where ( ) denotes the time average. That is 


(C.l) 


t2 


<>-r=r f ( ) 

^ 1 J tl 


dt 


Note also that 


£ _ P+P' _ ( p+p')(p-p') 

9 p+p' ( p+p')(p-p ') 

_ (p p - pp’ ± pp'- p’p') 

(p 2 - P' 2 ) 

= (p p - pp' + pp’- p>')(p 2 ±_/? ) 
(p 2 - p’ 2 ) (p 2 +p’ 2 ) 


(C.2) 
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The numerator becomes 

(p 2 -p' 2 ) (p 2 + p’ 2 ) = p 4 -p’ 4 
Neglecting p' 4 (small) we have 


2 = L 

' P 4 


pp 3 + (+ p 3 p' - pp 2 p’) + (- p 2 p'p' + ppp' 2 ) + (+ pp'p' 2 - pp’ 3 ) + 


Thus 


E = — -t £ + i v v- 
* r-l p5 Vj 


= e + e ( !> + e ( 2) + e ( 3 ) + e ( 4) 


where 


1 p i _ _ 

e = CTt + Z v j v J 

I n 


e ( i) = ^=r + v j v i 


p p* 


e < 2) = - x=r p'p’ - ^ ^’ 2 ) + 5 v j v j 


7 " A V ‘ • P 3 




7 ‘ i V P 4 

1 pV 3 

e (4 ) =~j=[ “ 

' P 

It follows from (C.2) through (C.6) that 

pE = (p + p') (e + e ( t) + e { 2) + e ( 3) + e ( 4) 

pE = pe -f t(pe ( 4 ) + p'e) + e 2 (pe ( 2) + p'e ( 4) ) 
+ f3 (P e ( 3) + P®( 2) ) + f 4 (P e ( 4) + P'®< 3 ) ) 

= 7=rP + j^i'i 


+ e 


+ e 2 


(p 1 - £ p’) + PVjVj + -^ £ p' + ^ V jV ]] 

1 p P J 

- A- (r p>' - !, p’ 2 ) + § Vi + i - L , p' 2 ) + p V 


7 — i * • -2 

1 P P 


7—1 -2 

' p p 


(- P> ,J ) 
(C.4) 


(C.5) 


(C.6) 
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= Gjlri-Vj) + e [t?t + 5 vi + f Vi 



(C.7) 


where the energy growth factor e was introduced with powers corresponding to the number 
of multiples of perturbed variables. 

Similarly, 

pVi E,i = (vj + v-) (p + p) (e + e ( l5 + e ( 2) + e, 3) + e ( 4) )j 

= Vipe,i ( c - 8 ) 


+ e pe ti v[ + vj (pe (1)>i + P%i) 

+ e 2 v[ (pe ( t) ,i + p’e.i) + (pe< 2) (i + p’e ( t) J 
+ e 3 v- (pe ( 2) ,i + p'e { 0 + Vj (pe ( 3 , ,i + p’e ( 2) .j) 

+ e 4 vj (pe { 3 ) ,i + p'e { 2) ti ) + v t (pe ( 4 , ^ + p'e ( 3) 

(C.9) 

S (pvj.i = (S ( t) + S ( 2) + S ( 3 ) + S ( 4) )(pvj + pv; + vjp' + pv;)^ 

= e (pvi),i S ( t ) 

+ f2 S( !) (pvjji + S( 2 ) (pvjji + S( J) (p Vj),ij 
+ f3 S( 3 ) (pVj),i + S, 2 ) (pvj.i + S( 4 ) (p ViJ.j + S( 2 ) (ViP ),i 
+ e 4 S (4) (pVi),i + S ( 3 ,(pvi) >i + S ( 3 ) (v i p , ), i + S (2 ) (p'v;),i 

(C.10) 


t~a 
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v iP-i = 


'£ + 

( E- 

-M 

1 + (-^- + 

V’) 

p 

p 

p 2 

p 2 



+ (^_L^) + (^) 


p 3 p 4 


▼iP.i + v iP’>i + V iP>i + V iP>i 


= E *iP»i 


+ € 


+ f 


^ (ViP'.i + v iP.i) + ViP’i ( E - \ p) 
p p p 

? (vip'.i) + ( £ - ^ P’) (ViP'.i + y'iPn) 
l P 


P P 2 


+ + *tP' 2 ) V>i 


P 2 P 3 


+ e* 


( e ’- E ;^) V ^i + (- e: 4 + E ;0 ,, )(V.. + »».J 

P P 

+ (^' ! + !V>V, 


p 2 p 3 


p 3 p 4 


+ f4 [(-^ + V' 2 ) v ^'-i + -^rp' 3 )(v , i + v ^.i) 

L p 2 p 3 p 3 p 4 

+ (_E^ 3 ) 

P 

PViVjVj,i = [p + p’][Vi + V-][Vj + vj][Vj,i + Vj, J 

= [pvj + p'vj + py[ + p’v^VjVj.i + vjVj.i + VjVj.i + vjvj^ 

= PViVjVj.i 

+ « [pVi( V jVj,i + VjV],j) + (pvj + pv-) VjVj.J 
+ e 2 [pViVjvj^ + (p'Vi + pv^vjvj.i + VjVj.i) + p'vJVjVj,^ 
+ e 3 [(p'Vi + pvj) VjVj.i + p'v; (VjVj.i + VjVj.i)] 

+ e 4 [p'v- vjVjjJ 


(C.ll) 


55 


p Vi E = pev i 

+ e [pev^ 4- vj (pe, t) + p'e)] 

+ e 2 [v](pe, d + p'e) + (pe, 2) + p'e, t ) )] 

+ e 3 [v](pe< 2) + p'e, 1} ) + Vj (pe, 3) + p'e, 2) )] 

+ 6 4 [v-(pe, 3) + p'e, 2) ) + vj (pe, 4) + p'e, 3) )] 


PVjS = € (pvjS, d ) 

+ e 2 [S, d pv] + S, j) pvi + S, i, p’vj 
+ e 3 [S, 3) pvi + S, 2 ,pv; + S, t) p'vi + S,2)ViP'] 

+ f 4 [S, 4) pVi + S,3, pVi + S, 2 )PVi+ S( 3 ) VjP ] 
p Vi = (p + p')(Vi + v])= pvj + e(pv- + p'Vi)+ e 2 p'vi 
a ij = ^ij + a ij 

where 

2 

^ij = - P^ij + M( v iij + v j>i) _ 5 M v k>lAj 
(?[■ = - p'^ij + m(v|, j + Vj,i) - 1 M v kik^ij 


Thus, 


a ij v j = (^ij + a ij)(*j + v j) = ^ij^j + £ (^ij v j + <jVj) + e Mj v j 


Substituting the above relations into (C.l) yields 
| f [ e ! E 1 + f ! E ! + e *E J ] = e 2 I 1 + e I I 2 + e , I 2 


where 


dQ> 


(C.12) 


(C. 13) 
(C.14) 
(C.15) 

(C. 16) 
(C.17) 

(C.18) 
• (C.lfl) 


E ' = <| n [i v i v i + ^i v i 


(C.20) 
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dn> 

<|^[v i (p e ( D ,i + p'e ti ) + Vi(p e ( 2) ,i + p'e { , £ ) 

f P 

+ S( D (pVj),i + S( 2) (pv-),i + S( D (p'vO.i — j 7 ( v i P'n) 

' n 


(C.21) 

(C.22) 


P P . . . , , p y p , 2 x - - 1 

P ) ViP.i- 

n J 


+ ( : P') (V»i + v iP>i) + (“ 2 + 

P p 2 P P 


- {p v i v j v j.i + (p’ v i + P v i)( v j v j>i + v j v j>i) + ^ v i v j v ju} dQ ) 


+ < [- {vi (pe ( D + p'e) + (pe ( 2) + p'e ( d )} 

-{S (1) pvi+ S ( 2 ) pvj} + {— p'^ij + P( v i>j + v],i) -§MVk, k ^j} • vj] nj dr) 

(C.23) 

(J o [{vi(pe ( 2) )i + p'e< 2) >i) + Vi (p e ( 3 ) ,i + pe { 2) ,i)} 

+ {S ( 3 ) (pvi),i + S ( 2) (P v i),i + SdjCp'v-J.i} 

P’ P P 'P P 

- {( : P') v iP’.i + (“ — +- P' 2 )(viP’.i + vip.i) 

P P 2 P 2 P 3 

P'P' 2 P _ _ I 

+ (— P ) ViP,i 

P 3 P 4 J 

- {(P’vi + pvi) vjvjji + p'vi (vjvj.i + VjVj.j] dfi > 

+ < | [- {vi (pe, 2) + p'e, u ) + Vi (pe ( „ + p'e ( 2) )} 


- {S ( 3 ) PVi + S ( 2 > pvi + S ( i, p’vi)] njdr) 


(C.24) 
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I 3 = ( { v i(P e ( 3 ) >i + P e ( 2 ) >i) + v i(P e ( 4) >i + p' e l 3) »i)} 


+ {S (4 , (pVi),i + S ( 3) (pvD.j + S ( 2 )(p’ v i))i} 


P P P 


PP^ P 


- {(— + — p' 2 ) vip'.i + ( — 7 P ,3 )(ViP',i + y'iPn) 

p 2 p 3 p 3 p 


+ (-~r - [p' v i v j v j’i]] dQ > 

p J 


+ ( [- { v i (p e ( 3) + P e ( 2 ) ) + V i (P e ( 4) + P e ( 3) }} 

J r 

- {S( 4) pVi + S( 3 , pv- + S( 2) p'vj}] riidr) 

Performing the differentiation as implied in (C.19), we obtain 

fa f21 1 + f3 ! 2 + f4l 3 


3t 


2fE t + 3e 2 E 2 -f* 4e 3 E 3 
E <> 2 E on 'i 


3E 

= ( eI i + f2l 2 + f3E 3) 5^ j 1 “ e 


+ t‘ 


9 2 . *^3 

Li ( ^-E7J 


(C.25) 


(C.26) 


where higher order terms and those terms much smaller than unity have been neglected. 
Thus, finally, we obtain 


de 

at 


- a { e - a 2 e 2 - a 3 e 3 = 0 


(C.27) 
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APPENDIX D 

INTEGRANDS OF E t , E 2> E 3 , I t , I 2 , 1 3 


a ( i> =f vjvj + p’vjv: 


a< 2) =^' v j v j 


a ( 3) = £-p ' 4 
P 4 


b ( 1> = (v£ + pv i )(-i r (§- ^ p') + v j v]) Ji + v-p' + (i ? + 5 VjVj),i 


^7— I -2 

P p 


+ pv i (- pr ^ ^ + 5 v J v j)’i + (K + (^pr ?- t=t £ ) 

1 p z p 6 1 ~ 1 


7-1 - 7- _ 

' P P 


+ (P *i).i (- j) (A- ( E ’) 2 - A ( £ Y) - E Mp.;) - ( E - ^ p')(v,ft ; 

' P 1 P P P P 


+ v[p, i) - (- ^ + ^7 P' 2 ) ViP.i - ^i v j v j>i - (p'v i + P v i) ( v j V + VjVj.i) 


p 2 p 3 


pVj VjVj.i 


c“> =(vip + p l v i )(-l T (^-£-p l ) + v j v]) + vip' (-Ij- ^ + jVjVj) 


^ 7— i ' _ - 2 

' P P 


+ PVi P' 2 ) + 5 V j V j)u + (K + P'Vi),! (A ? 


7 " A V P 3 


+ (p y i) (-5) (p4 (^) 2 - ^r(^) 2 ) + p'^ij — P ( v i’j + V j>i) 

’ - ’ p > p 


+ f M v k>k v j] n j 


A 
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b< 2) = (v[p + p’Vi)(-i ~\P 2 ) + 5 V j V j)’i + P ' V i (^=T ^ 

1 p 2 r p p 

+ v j vj), i + pv i (^ r (^ — ?7 P' 3 ))>i + ^ ^=T (?") (M’i 


v 7 -i -3 -4 

1 P P 


A( e V-i( £ ') 2 ] (p v; + pvj.i + 


I 7-i ' 7 

L ' P P 


7-1 - 7~1 ^ 

1 P ' P 


(P' v i).i 


' .’2 


( £ “ ^ p) y\p\ i - (- ^ + *7 p' 2 ) ( v iP'n + v'iP’i) - - ?7 P*) v i^’i 


P P 2 


P 2 P 3 


P 3 P 4 


(p'Vi + pv[) vjvj fi - p'v; (vjVj.i + VjVj.i) 


c< 2) = (v|p + p’Vi) (- -^j- p' 3 ) + \ VjVj) + p'v- (i fi- - ^ P 1 ) 


7 -i -2 -3 

' P P 


P P 


+ VjVj) + pv i (^i - ?7 P' 3 )) + l lyir " J=1 ^ 


^ V P 4 


1 


(jf ) 2 - £ Cr) 7 ] (p y i + p' y i) + [ T ^; jr - p p’ v a “i 


b 131 =(»;? + »/) ^ (^ -^p' ! ),i + p'v;((-^i)(^p> - -^p' ! ) 


y 


+ J Vjyj),, - p V, ^ J [^r (^) 3 - ^ t^) 3 ] (pv! + P'Vi),, 


- 2 & (f ) 2 - & (f) 2 ) (p'vi), -15 l^r (fy - £ (p‘) (pv,),. 


l 7 -l ' 7 _ 

1 P P 


> ' /i 1 2 sv3 


_ ( _t£ + E-p' 2 ) vip',, - (E-£ -££ KviP',, + v;p,,)-(-£E )»,?„ 
p 2 P 3 P 3 P 4 P 


P V i V j V j’i 
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>( 3 ) 


= (v* p + v 1 ) i ~^~a p '^ + pv i ((“ t=t) (~ 2 vp - ^3 p 
7 9 l P p P 


+ 5 v j v j) - PVi ^ ( £ ^' 3 ) + J [^r ( £ ') 3 - (f) 3 ](P v i + P'^i)’i 


p 

1 r 1 


i (i ( E ) 2 - ? ) 2 ) - (p-i) - k ( E ) 4 - ( E ) 4 1 ^ n i 


APPENDIX E 


Listing of Computer Program (ECI-1) 


PROGRAM TG1D 

C 

PARAMETER (NELEM=200 , NPOIN=201 ) 

CALL DINPUT 
CALL LPMASS 
CALL ITERAT 

C 

STOP 

END 

C 

SUBROUTINE DINPUT 

PARAMETER (NELEM=200 , NP0IN=201 ) 

PARAMETER ( NNODP=2 , NGAUS=2 , NCONS=3 ) 

COMMON/AAAA/A(NPOIN) ,DA(NPOIN) 

COMMON/DOMA/DO , UO , EO , PO 

COMMON/INIT/DENSY (NPOIN) .UVELY(NPOIN) , ENEGY (NPOIN) , PRESY(NPOIN) 
COMMON/COOR/XX (NPOIN) , LNODS (NELEM , NNODP ) 

COMMON/PRTY/CAPAV , CAPAP , CGAM , CONDT . VISCY 
COMdON/TIME/CFLNB , DTI ME , ITMAX 
DIMENSION XI(IOO) , AI(100) 

DIMENSION EA( NELEM) , EDA (NELEM) 

C 

C READ IN' FLOW PROPERTIES AND TEMPORAL PARAMETERS 
C 

REaO (19,*) ITMAX 
IRE AD=2 
CGAM=1 . 22 

CArw=i.o 

C0’OT=0.0 

vis:y=o.o 

CFLNB=0. 6 
DTI\1E=0 . 025 

CC WRITE (6, 2010) CGAM , CAPAV , CONDT , VISCY 

C WRITE (6, 2020) CFLNB . DTI ME , ITMAX 

C 

C READ IN NODAL CONNECTIVITIES 
C 

C WRITE (6, 2030) 

DO 10 1=1, NELEM 
LNODS ( I , 1 ) =1 
LNODS ( I , 2 ) =1+1 

C WRITE(6, 1000) I , LNODS (1,1) , LNODS (1,2) 

10 CONTINUE 

1000 FGRvtAT(lX, 15 , 5X, 215) 

C 

C READ IN NODAL COORDINATES 
C 

CCRD=0. 0254*5 . 1527 
XTH=0 . 0 

ATH=3 . 14*C0RD* *2 
IN-82 

DO 151 1=1, IN 

READ (17,152) II,XI(I) ,AI(I) 

C PRINT*, I, XI (I) ,AI(I) 

XI ( I ) =CORD*XI ( I ) 

AI < I ) =ATH*AI ( I ) * * 2 
C PRINT*, I,XI(I) , AI ( I ) 

151 CONTINUE 

152 FORMAT (5X,I5,2X,2E12.5) 


DX= (XI(IN)-XI(l) )/ FLOAT (NELEM) 

XX(1)=XI (1) 

A ( 1 ) =AI ( 1 ) 

XX(NPOIN)=XI (IN) 

A(NPOIN)=AI (IN) 

DO 20 1=1 , NPOIN-1 
XX ( I +1 ) =XX ( I ) +DX 
20 CONTINUE 

DO 153 1=2, NPOIN-1 
XA=XX ( I ) 

DO 154 J=1 , IN-1 

IF (XA . GE . XI ( J ) . AND. XA. LE. XI ( J-t-1) ) THEN 
SLOPE=(AI ( J+l) -AI ( J) ) / (XI ( J+1)-XI ( J) ) 

A ( I ) =AI ( J ) +SLOPE* ( XA-XI ( J ) ) 

END IF 

154 CONTINUE 
153 CONTINUE 

DO 5003 1=1, NELEM 

EA ( I ) =0 . 5 * ( A ( I ) +A ( 1 + 1 ) ) 

DXX=XX ( 1+1 ) -XX ( I ) 

EDA ( I ) = ( A ( 1+1 ) -A ( I ) )/DXX 

5003 CONTINUE 

DO 5004 1=2, NPOIN-1 

DA ( I ) =0 . 5* ( EDA ( 1-1 ) +EDA ( I ) ) 

5004 CONTINUE 

DA ( 1 ) =EDA ( 1 ) 

DA (NPOIN)=EDA( NELEM) 

DO 85 1=1 , NPOIN 

C WRITE (6, 2500) I , XX ( I ) , A ( I V , DA ( I ) 

85 CONTINUE 

2500 FORMAT (1X,I5,F10.5,2E15.5) 

C 

C READ IN INITIAL CONDITIONS 
C 

C AIR PROPERTIES @ T=1000 K 

REC=2 . 67E+5 
CGAM=1 . 2 
CAPAV=1 . 7 
C CAPAV=1000 . 

READ (19,*) APRES , ATEMP 
PATM=AP RE S / 1 4 . 7 
PSTG=PATM*9 . 8E4 
CTEM= ( ATEMP-460 . )*5./9.+273 
C CTEM=1000 . 

CMACH=0 . 2 

CGAS=1 . 987*1000 .*4.184 
CGRAV=9 . 8 

CWGT=0. 79*28 . +0.21*32. 

CSND=SQRT ( CGAM*CTEM * CGAS /CWGT ) 

' CVEL=CMACH*CSND 
CSQR=0 . 5*CVEL**2 
CAPAV=CAPAV*CGAS/CWGT 
CENG=CAPAV*CTEM+CSQR 
CPRE=PSTG 

CRHO=CPRE/ ( ( CGAM-1 . ) * (CENG- CSQR) ) 
CENT=CENG+CPRE/CRHO 

C PRINT* .CSND.CVEL, CRHO, CENG, CPRE.CAPAV 

C £TOP 

C 

CAPAP=CGAM*CAPAV 
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DO=CRHO 

PO=CPRE 

UO=CVEL 

V0=0 . 0 

TO=CTEM 

EO=CENG 

HO=CENT 

CAPAP=CGAM*CAPAV 
DO 30 1=1 , NPOIN 
PRESY ( I ) =P0 
UVELY ( I ) =0 . 0 

ENEGY ( I ) =CAPAV*T0+0 . 5*UVELY( I ) **2 

DENSY ( I ) =PRESY ( I ) / ( (CGAM-1 .0) * (ENEGY ( I ) -0 . 5*UVELY ( I ) **2) ) 
30 CONTINUE 
PRESY ( 1 ) =P0 
UVELY ( 1 ) =U0 
ENEGY ( I ) =E0 
DENSY ( 1 ) =D0 
C 

C RESTART PROCEDURES 
C 

IF ( I READ . EQ . 1 ) THEN 

READ(11, 1060) <XX( I ) , 1=1 , NPOIN) 

READ( 11 , 1060) ( A( I ) , 1=1 , NPOIN) 

READ (11 , 1060) ( DENSY ( I ) , 1=1 .NPOIN) 

READf 11 , 1060) (UVELY( I ) , 1=1 , NPOIN) 

READ (11 , 1060) ( ENEGY ( I ) , 1=1 , NPOIN) 

READ (11 , 1060) ( PRESY ( I ) , 1=1 . NPOIN) 

END IF 

1060 FORMAT (5 ( 200 ( 4E15 . 5 , / ) ) ) 

C 

C WRITE OUT COORDINATES AND INITIAL CONDITIONS 

C 

C 

RETURN 

C 

2010 FORMAT (//' PHYSICAL PROPERTY' ,/ ' *****************' // 

' CGAM = ' , F7 . 4 , 4X , CAPAV *',F7.4,4X, 

' CONDT =' . F7 . 4 , 4X , VISCY =',F7.4) 

2020 FORMAT*//' INITIAL TIME STEP',/' * ^ *************** ' ,// 

' CFLNB =' , F7 . 4 , 4X , ' DTI ME =',F7.4,4X, 

' ITMAX =' , 15) 

2030 FORMAT (//' ELEMENT TOPOLOGY',/' ****************' ,// 

' ELEMENT' , 6X, 'NODE NUMBERS') 

2040 FORMAT*//' NODE POINT DATA',/' ***************' ,// 

, 5X , 'NODE' , IX, 'X' , 10X, 'DENSY' ,5X, 

'UVELY' , 5X, 'ENEGY' , 5X, ' PRESY' ) 

2045 FORMAT (5X, I4,3X,5(F7,4, 3X) ) 

C 

END 

C 

SUBROUTINE LPMASS 
PARAMETER ( NELEM=200 ,NPOIN=201) 

PARAMETER (NNODP=2) 

COMMON/COOR/'XX ( NPOIN) . LNODS ( NELEM , NNODP ) 
COMMON/MASS/GMASS ( NPOIN ) 

DIMENSION F I ( 2 ) , POSGP ( 2 ) , WE IGP ( 2 ) 

INITIALIZATION OF LUMPED MASS 



DO 10 1=1 , NPOIN 
GMASS ( I ) =0 . 0 
10 CONTINUE 
C 

C*** SET UP POSITIONS AND WEIGHTS FOR 2 POINT GAUSS RULE 

C 

P0SGP(1)=0. 5773502691 
POSGP ( 2 ) =-POSGP ( 1 ) 

WEIGP ( 1 ) =1 . 0000000000 
WEIGP ( 2 ) =WEIGP f 1 ) 

C 

C ASSEMBLE LUMPED MASS 
C 

DO 100 IELEM=1 , NELEM 

SLETH=ABS ( XX ( LNODS ( IELEM , 2 ) ) -XX ( LNODS ( IELEM , 1 ) ) ) 

C 

C INTEGRATIONS 

C 

DO 90 IGAUS=1 , 2 
C 

D JA=0 . 50* SLETH* WEIGP ( IGAUS ) 

XI=POSGP ( IGAUS ) 

FI(1)=0.50*(1.0-XI) 

FI ( 2 ) =0 . 50* ( 1 . 0-t-XI ) 

C 

DO 30 1=1 , NNODP 
K=LNODS( IELEM, I) 

SHAPX=FI ( I ) 

DO 30 J=l, NNODP 
SHAPY=FI ( J ) 

GMASS (K) =GMASS ( <) +SHAPX*SHAPY*DJA 
30 CONTINUE 
90 CONTINUE 
100 CONTINUE 
C 

C STORE IN THE OUTER- CORE MEMORY 
C 

RETURN 

END 

C 

SUBROUTINE SDTIME 

PARAMETER ( NELEM=200 , NP0IN=201 ) 

PARAMETER (NNODP=2) 

COMMON/ AAAA/ A < NPOIN) ,DA(NPOIN) 

COMMON/AREA/ AREAL (NELEM) 

COMMON/COOR/XX ( NPOIN ) , LNODS ( NELEM, NNODP ) 

COMMON / INIT/DENSY ( NPOIN ) , UVELY ( NPOIN ) , ENEGY ( NPOIN) , PRESY(NPOIN ) 
COMMON/ PRTY/CAPAV , CAPAP , CGAM , CONDT , VISCY 
COMMON/TIME/CFLNB , DTIME , ITMAX 
DIMENSION TIMEL( NELEM) 

DIMENSION DENSM( NNODP) , PRESM ( NNODP ) , UVELM ( NNODP ) , WELM ( NNODP ) 

C 

C EVALUATE TIME STEP IN EACH ELEMENT 
C 

DO 10 IELEM=1, NELEM 
DO 20 J=l, NNODP 
K=LNODS< IELEM, J) 

DENSM ( J ) =DENSY ( K) 

UVELM ( J ) =UVELY ( K ) 

PRESM ( J)=PRESY(SC) 
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20 CONTINUE 
DABSM=0 . 0 
UABSM=0 . 0 
PABSM=0 . 0 
AA=0 . 0 

DO 30 1=1 , NNODP 
DABSM=DABSM+0 . 5 *DENSM ( I ) 

UABSM=UABSM+0 . 5 *UVELM ( I ) 

PABSM=PABSM+0 . 5*PRESM( I ) 

AA=AA+0 . 5 * A ( LNODS ( IELEM , I ) ) 

30 CONTINUE 
C 

SLETH=ABS ( XX ( LNODS ( I ELEM , 2 ) ) -XX ( LNODS ( IELEM , 1 ) ) ) 

UVABS=ABS (UABSM) 

CSPED=SQRT (CGAM*ABS ( PABSM) /ABS(DABSM) ) 

TIMEL ( IELEM) =CFLNB * SLETH/ (UVABS+CSPED) 

C TIMEL ( IELEM) =CFLNB*SQRT (AREAL ( IELEM) ) / (UVABS+CSPED) 

10 CONTINUE 

C 

C FIND MINIMUM TIME STEP 

C 

DTIME=TIMEL ( 1 ) 

C CFLLL=TIMEL ( 1 ) 

DO 40 IELEM=2 , NELEM 

IF (TIMEL (IELEM) .LT.DTIME) DTIME=TIMEL ( IELEM) 

C IF(TIMEL ( IELEM) .GT.CFLLL) CFLLL=TIMEL ( IELEM) 

40 CONTINUE 

PRINT* , 'CFLNUMBER === , CFLLL 

RETURN 
END 

C 

SUBROUTINE MATRIX ( I ITER , IEQNS , IELEM) 

PARAMETER (NELEM=200 , NPOIN=20I ) 

PARAMETER ( NN0DP=2 , NEQNS=3 , NGAUS=2 ) 

COMMON/AAAA/A(NPOIN) ,DA(NPOIN) 

COMMON/DOMA/DO , UO , EO , PO 
C0MM0N/BCBC/D1 , U1 , El , PI 
COMMON/AREA/AREAL (NELEM) 

COMMON /COOR/XX ( NPOIN ) , LNODS ( NELEM , NNODP ) 

COMMON/ PRTY/CAPAV , CAPAP , CGAM, CONDT , VISCY 
COMMON/TIME/CFLNB , DTI ME , ITMAX 

COMMON/ INI T/DENSY ( NPOIN ) , UVELY( NPOIN ) , ELEGY (NPOIN) , PRESY ( NPOIN ) 
COMMON/HALF/DENSH( NELEM) , UVELH ( NELEM ) , ENEGH ( NELEM) , PRESH ( NELEM) 
COMMON / EQNS / EQRHR (NPOIN) , EQRHU(NPOIN) , ECRHE (NPOIN) 

DIMENSION POSGP ( NGAUS ) , WEIGP ( NGAUS ) ,FI<2) ,DX(2) 

DIMENSION UHALF (NEQNS) , FHALF ( NEQNS ) , FLUXH ( NEQNS ) , RHALF ( NEQNS ) 

C 

c *** SEX UP POSITIONS AND WEIGHTS FOR 2 POINT GAUSS RULE 
C 

AH=0 . 5* ( A ( LNODS ( IELEM , 1 ) ) +A ( LNODS ( IELEM , 2 ) ) ) 

DAH=0 . 5 * ( DA ( LNODS ( I ELEM . 1 ) ) +DA ( LNODS ( I ELEM , 2 ) ) ) 

POSGP(1)=0. 5773502691 
POSGP ( 2 ) =-POSGP ( 1 ) 

WEIGP (1)=1 .0000000000 
WEIGP ( 2 ) =WEIGP ( 1 ) 

C 

C LOOP TO CARRY OUT GAUSS INTEGRATION 

C NOTE : PERFORMED JUST ONCE IN EACH TEMPORAL ITERATION 
C 
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IF ( IEQNS . NE . 1 ) GO TO 20 
C 

AREAL ( IELEM) =ABS (XX (LNODS ( IELEM , 2 ) ) -XX ( LNODS ( IELEM , 1 ) ) ) 

DO 10 J=1 , NEQNS 
RHALF ( J ) =0 . 0 
10 CONTINUE 

EVALUATE INTEGRATIONS IN THE RIGHT-HAND SIDE OF HALF STEP 

DO 70 IGAUS=1 , NGAUS 

DJA=0. 50* AREAL (IELEM) *WEIGP ( IGAUS ) 

DTA=0. 50*AREAL( IELEM) 

XI=POSGP( IGAUS) 

FI (1)=0. 50* (1 . 0-XI) 

FI (2 ) =0 . 50* ( 1 . 0+XI ) 

DX ( 1 ) =-0 . 50/DTA 
DX ( 2 ) = 0. 50/DTA 

EVALUA'I E PREVIOUS VARIABLES AND FLUXES AT GAUSS POINTS 
IN THE HALF STEP 

S0RCE=0 . 0 
DO 40 J=l, NEQNS 
UHALF ( J ) =0 . 0 
FH/*LF ( J ) =0 . 0 
40 CONTINUE 

DO 50 1=1 , NNODP 
K=l NODS ( IELEM, I) 

UHALF ( 1 ) =UHALF ( 1 ) +DENSY (K) * A(K) *FI ( I ) 

UHALF ( 2 ) =UHALF ( 2 ) +DENSY (K) *UVELY (K) *A(K) *FI( I) 

UHALF ( 3 ) =UHALF ( 3 ) +DENSY ( K) *ENEGY ( K) *A ( K) *FI ( I ) 

FHALF ( 1 ) =FHALF ( 1 ) +DENSY ( K) *UVELY(K) *A(K) *DX(I) 

FHALF ( 2 ) =FHALF ( 2 ) + (DENSY ( K) *UVELY ( K) * *2+PRESY (K) ) *A(K) *DX ( I ) 

FH, LF < 3 ) =FHALF ( 3 ) +UVELY (K.) *A (K) * (DENSY ( K) *ENEGY (K) +PRESY ( K) ) 

*DX ( I ) 

SORCE=SORCE+PRESY (K)*DA(K)*FI(I) 

50 CONTINUE 

ORGANIZE RIGHT-HAND SIDE OF HALF STEP 

RKALF ( 1 ) =RHALF ( 1 ) +DJA* (UHALF ( 1 ) -0 . 5 *DTIME*FHALF ( 1 ) ) 

RHaLF ( 2 ) = RHALF ( 2 ) +DJA* (UHALF ( 2 ) +0 . 5 *DTIME* ( SORCE-FHALF ( 2 ) ) ) 

RHALF ( 3 ) =RHALF ( 3 ) +DJA* (UHALF (3 ) -0 . 5 *DTIME* FHALF ( 3 ) ) 

70 CONTINUE 
C 

C CALCUI ATE EACH VARIABLE AT THE HALF STEP 
C - . 

^ DENSH( IELEM) =RHALF ( 1 )/( AREAL ( IELEM) *AH ) 

UVELH ( IELEM) =RHALF (2 ) / (DENSH ( IELEM) * AREAL ( IELEM) *AH) 

ENEGH ( IELEM) =RHALF ( 3 ) / ( DENSH ( IELEM) * AREAL ( IELEM) *AH ) 

PRESH ( IELEM) =(CGAM-1 . 0 ) *DENSH ( IELEM) 

* (ENEGH ( IELEM ) -0 . 5* UVELH ( IELEM) *UVELH ( IELEM) ) 

20 CONTINUE 

C 

C CALCULATE FLUX TERMS AT THE HALF STEP 

C a 

FLUXH(1)=DENSH( IELEM) *UVELH ( IELEM ) *AH 

FLUXH ( 2 ) = ( DENSH ( IELEM ) *UVELH ( IELEM) *UVELH ( IELEM ) -i-PRESH ( IELEM ) > * AH 
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FLUXH ( 3 ) =UVELH ( IELEM ) * ( DENSH ( I ELEM ) * ENEGH ( IELEM ) +PRESH ( IELEM ) ) * AH 
SORCH =PRESH ( IELEM) *DAH 

EVALUATE INTEGRATIONS IN THE RIGHT-HAND SIDE AT THE FULL STEP 
DO 80 IGAUS=1 , NGAUS 

DJA=0. 50*AREAL( IELEM) *WEIGP < IGAUS) 

DTA=0 .50* AREAL (IELEM) 

XI=POSGP( IGAUS) 

FI (1)=0. 50* (1 . O-XI) 

FI (2) =0. 50* (1 . O+XI) 

DX ( 1 ) =-0 . 50/DTA 
DX ( 2 ) = 0. 50/DTA 
C 

C EVALUATE RIGHT-HAND SIDE AT THE FULL STEP 

C 

DO 110 1=1 , NNODP 
K=LNODS( IELEM, I) 

CARXI=DX ( I ) *DJA*DTIME 
C GOTO (111,112,113), IEQNS 

111 EQRHR ( K ) =EQRHR ( K ) +FLUXH ( 1 ) *CARXI 

C GO TO 110 

112 EQRHU(K) =EQRHU(K) +FLUXH (2) *CARXI+SORCH*FI ( I ) *DJA*DTIME 

C GO TO 110 

113 EQRHE ( K) =EQRHE ( K ) +FLUXH ( 3 ) *C ARX I 
110 CONTINUE 

80 CONTINUE 

C 

RETURN 

END 

C 

SUBROUTINE BDFLUX( IEQNS) 

PARAMETER ( NELEM=200 , NPOIN=201 , NN0DP=2 ) 

COMMON/AAAA/A ( NPOIN ) ,DA(NPOIN) 

COMMON /COOR/ XX ( NPOIN ) , LNODS ( NELEM , NNODP ) 

COMMON / PRTY/ CAP AV , CAPAP , CGAM , CONDT , VISCY 
COMMON/TIME/CFLNB , DTI ME , ITMAX 

COMMON/ INIT/DENSY( NPOIN) , UVELY ( NPOIN ) , ENEGY (NPOIN) , PRESY(NPOIN) 
COMMON/HALF/DENSH( NELEM) , UVELH ( NELEM) , ENEGH ( NELEM) , PRESH (NELEM ) 
COMMON/EQNS/ EQRHR (NPOIN) , EQRHU( NPOIN ) , EQRHE ( NPOIN) 

C 

C EVALUATE AVERAGES INSIDE DOMAIN ELEMENT 

C 

K11=LN0DS (1,1) 

K12=LN0DS (1,2) 

KN1=LN0DS ( NELEM , 1 ) 

KN2=LN0DS ( NELEM , 2 ) 

C 

DRDA1=0 . 5 *DENSY (Kll) * UVELY (Kll ) * A ( Kll ) 

+0. 5*DENSY(K12) *UVELY(K12) *A(K12) 

DUDA1 =0 . 5 * ( DEN SY < K1 1 ) * UVELY ( K1 1 ) * UVELY ( K1 1 ) +PRESY ( K1 1 ) ) * A ( K1 1 ) 

+0. 5 * ( DENSY ( K12 ) * UVELY ( K12 ) *UVELY( K12 ) +PRESY ( K12 ) ) *A(K12) 
DEDA1=0 . 5 * (DENSY ( Kll ) * ENEGY ( Kll ) +PRESY ( Kll ) ) *UVELY ( Kll ) * A ( Kll ) 
+0.5* (DENSY ( K12 ) *ENEGY ( K12 ) +PRESY (K12 ) ) *UVELY ( K12 ) *A ( K12 ) 

C 

DRDAN=0 . 5 * DENSY (KN1 ) *UVELY ( KN1 ) *A( KN1 ) 

+0 . 5 * DENSY ( KN2 ) * UVELY ( KN2 ) *A ( KN2 ) 

DUDAN=0 . 5 * ( DENSY ( KN1 ) *UVELY i KN1 ) *UVELY ( KN1 ) +PRESY ( KN1 ) ) * A ( KN1 ) 

+0 . 5 * ( DENSY ( KN2 ) *UVELY ( KN2 ) *UVELY ( KN2 ) +PRESY ( KN2 ) ) *A ( KN2 ) 
DEDAN=0 . 5 * ( DENSY ( KN1 ) * ENEGY ( KN1 ) +PRESY ( KN1 ) ) *UVELY (KN1) *A(KN1 ) 
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+0.5* ( DENSY (KN2 ) *ENEGY ( KN2 ) +PRESY ( KN2 ) ) *UVELY(KN2) *A(KN2 > 


EVALUATE BOUNDARY TERMS AT THE HALF STEP 


AH1=0 . 5 * ( A ( Kll ) +A ( K12 ) ) 

AHN=0 . 5* ( A ( KN1 ) +A ( KN2 ) ) 

DRDH1=DENSH ( 1 ) *UVELH ( 1 ) *AH1 

DUDH1= ( DENSH ( 1 ) *UVELH ( 1 ) *UVELH ( 1 ) +PRESH ( 1 ) ) *AH1 
DEDH1= ( DENSH ( 1 ) *ENEGH ( 1 ) +PRESH ( 1 ) ) *UVELH ( 1 ) *AH1 

C 

DRDHN=DENSH ( NELEM ) *UVELH (NELEM) *AHN 

DUDHN= (DENSH (NELEM) *UVELH (NELEM) *UVELH ( NELEM) +PRESH < NELEM) ) * AHN 
DEDHN=( DENSH (NELEM) *ENEGH ( NELEM) +PRESH ( NELEM) ) *UVELH ( NELEM) *AHN 
C 

C ZERO-TH TIME STEP 
C 

DRDN1=DENSY ( 1 ) *UVELY ( 1 ) * A ( 1 ) 

DUDN1= ( DENSY ( 1 ) *UVELY( 1) *UVELY ( 1 ) +PRESY ( 1 ) ) *A( 1 ) 

DEDN1= (DENSY ( 1 ) *ENEGY ( 1 ) +PRESY ( 1 ) ) *UVELY ( 1 ) * A ( 1 ) 

C 

DRDNN=DENSY (NPOIN) *UVELY ( NPOIN) *A(NPOIN) 

DUDNN= ( DENSY ( NPOIN ) *UVELY ( NPOIN) *UVELY ( NPOIN ) +PRESY ( NPOIN ) ) 

* A ( NPOIN ) 

DEDNN=(DENSY (NPOIN) *ENEGY ( NPOIN ) +PRESY ( NPOIN ) ) *UVELY ( NPOIN ) 
*A(NPOIN) 


C 

C 

C 

C 


C 


C 


C 


INCLUDE BOUNDARY GRADIENT TERMS INTO RHS VETCOR 


31 


32 


33 

30 


GO TO (31,32,33) , IEQNS 

EQRHR(l) =EQRHR( 1) -DTIME* ( -DRDN1-DRDH1 +DRDA1 ) 
EQRHR(NPOIN) =EQRHR ( NPOIN ) +DTIME* ( -DRDNN-DRDHN+DRDAN) 
GO TO 30 

EQRHU ( 1 ) =EQRHU ( 1 ) -DTIME* ( -DUDN1-DUDH1+DUDA1 ) 

EQRHU ( NPOIN ) =EQRHU (NPOIN) +DTIME* ( -DUDNN-DUDHN+DUDAN) 
GO TO 30 

EQRHE ( 1 ) ==EQRHE ( 1 ) -DTIME* ( -DEDN1 -DEDH1 +DEDA1 ) 

EQRHE( NPOIN) =EQRHE( NPOIN ) +DTIME* ( -DEDNN-DEDHN+DEDAN) 
CONTINUE 


RETURN 

END 


C 

SUBROUTINE SOLVER ( IEQNS , DELTA) 

PARAMETER ( NELEM=200 , NPOIN=201 ) 

PARAMETER ( NNODP=2 , NCONS=3 ) 

COMMON/MASS/GMASS ( NPOIN) 

COMMON/ AREA/ AREAL (NELEM) 

COMMON/COOR/XX ( NPOIN ) , LNODS ( NELEM , NNODP ) 

COMMON / EQN S /EQRHR (NPOIN) , EQRHU (NPOIN) , EQRHE ( NPOIN) 

DIMENSION DELTA (NPOIN) , EQRHS (NPOIN) , CDUMY (NPOIN) , GDUMY ( NPOIN ) 
DIMENSION POSGP ( 2 ) , WEIGP ( 2 ) , FI ( 2 ) 

C 

C*** SET UP POSITIONS AND WEIGHTS FOR 2 POINT GAUSS RULE 
C 

POSGP ( 1 ) =0 . 5773502691 
POSGP ( 2 ) =-POSGP ( 1 ) 

WEIGP ( 1 ) =1 . 0000000000 
WEIGP ( 2 ) =WEIGP ( 1 ) 

C 

GO TO (1,2,3) , IEQNS 
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1 DO 5 1=1 , NPOIN 

5 EQRHS ( I ) =EQRHR ( I ) 

GO TO 9 

2 DO 6 1=1, NPOIN 

6 EQRHS ( I ) =EQRHU ( I ) 

GO TO 9 

3 DO 7 1=1, NPOIN 

7 EQRHS ( I ) =EQRHE ( I ) 

9 CONTINUE 

C 

C READ LUMPED MASS FROM STORED TAPE 

C 

C SOLUTION PROCEDURE OF ALGEBRAIC EQUATIONS USING EXPLICIT 
C METHOD 
C 

C - LUMPED MASS 
C 

IF (NCONS . EQ . 1 ) THEN 
C 

DO 200 1=1, NPOIN 

DELTA ( I ) =EQRHS ( I ) /GMASS ( I ) 

200 CONTINUE 
END IF 

C 

C - JACOBI ITERATIONS 

C 

IF ( NCONS . EQ . 3 ) THEN 
DO 100 IC0NS=1. NCONS 
IF ( ICONS. NE. 1) GO TO 20 
DO 10 1=1, NPOIN 
10 GDUMY ( I ) =0 . 0 
20 CONTINUE 

DO 30 1=1, NPOIN 
30 CDUMY ( I ) =0 . 0 
C 

C COMPUTATION OF M*DU 

C 

DO 80 IELEM=1 , NELEM 

C 

C LOOP TO CARRY OUT GAUSS INTEGRATION 
C 

DO 70 IGAUS=1 , 2 
C 

DJA=0 . 50* AREAL ( 1ELEM) *WEIGP ( IGAUS ) 

DTA=0 . 50*AREAL ( IELEM) 

XI=POSGP( IGAUS) 

FI ( 1 ) =0 . 50* ( 1 . 0 XI) 

FI (2)=0. 50* (1 . O+XI ) 

C 

GINTP=0 . 0 
DO 50 1=1 , NNODP 
K=LNODS< IELEM, I ) 

GINTP=GINTP+GDUMY ( K) *FI ( I) 

50 CONTINUE 

DO 60 1=1, NNODP 
K=LNODS( IELEM, I ) 

CDUMY (K)=CDUMY(K.) +GINTP*FI (I ) *DJA 
60 CONTINUE 
70 CONTINUE 
80 CONTINUE 
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CALCULATION OF DELTA IN EVERY ITERATION 
DO 90 1=1 , NPOIN 

DELTA ( I ) = ( EQRHS ( I ) -CDUMY ( I ) > /GMASS ( I ) +GDUMY ( I ) 
90 CONTINUE 

DO 110 1=1, NPOIN 
110 GDUMY ( I ) =DELTA ( I ) 

100 CONTINUE 
END IF 

RETURN 

END 


SUBROUTINE LAPDUS ( IEQNS) 

PARAMETER ( NELEM=200 , NP0IN=201 ) 

PARAMETER (NNODP=2) 

COMMON /AAAA/ A ( NPOIN) ,DA(NPOIN) 

COMMON/ AREA/ AREAL ( NELEM ) 

COMMON/ I NIT/DENSY ( NPOIN ) , UVELY ( NPOIN ) , ENEGY ( NPOIN ) , PRESY (NPOIN) 
COMMON/EQNS/EQRHR( NPOIN) , EQRHU (NPOIN ) , EQRHE ( NPOIN) 
COMMON/COOR/XX (NPOIN) , LNODS ( NELEM , NNODP ) 

COMMON/TIME/CFLNB , DTI ME , ITMAX 

DIMENSION X(2) ,U(2) ,FI(2) ,DX(2) , POSGP ( 2 ) WEIGP (2) 

C 

C*** SET UP POSITIONS AND WEIGHTS FOR 2 POINT GAUSS RULE 

C 

POSGP(l) =0. 5773502691 
POSGP ( 2 ) =-POSGP ( 1 ) 

WEIGP ( 1 ) =1 . 0000000000 
WEIGP ( 2 ) =WEIGP ( 1 ) 

C 

C COMPUTATION OF ARTIFICIAL VISCOSITIES USING '.APIDUS' CONCEPT 

C 

DO 100 IELEM=1, NELEM 

C 

C ARTIFICIAL VISCOSITIES 

C 

DO 10 1=1, NNODP 
K=LNODS ( IELEM , I) 

X ( I ) =XX ( K.) 

U ( I ) =UVELY ( K) 

10 CONTINUE 
C 

DUDXA=ABS ( ( U ( 2 ) — U ( 1 ) ) / ( X ( 2 ) -X ( 1 ) ) ) 

C • . 

C 4.00P TO CARRY OUT GAUSS INTEGRATION 
C 

DO 100 IGAUS=1 , 2 

C 

DJA=0. 50*AREAL( IELEM) *WEIGP ( IGAUS ) 

DTA=0 . 50*AREAL ( IELEM) 

XI=POSGP( IGAUS) 

;.pi(i)=o.50*(i.o-xi) 

HI ( 2 ) =0 . 50* ( 1 . O+XI ) 

DX ( 1 ) =-0 . 50/DTA 
DX ( 2 ) = 0. 50/DTA 
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ddrdx=o . o 

DRUDX=0 . 0 
DREDX=0 . 0 
DO 40 1=1 , NNODP 
K=LN0DS ( IELEM , I ) 

DDRDX=DDRDX+DEN SY ( K. ) *A(K) *DX(I) 
DRUDX=DRUDX+DENSY (K) *UVELY(K) *A(K) *DX(I) 
DREDX=DREDX+DENSY(K) *ENEGY ( K) *A(K) *DX(I) 

40 CONTINUE 

ARTIFICAL VISCOSITY 

C0NSX= 1 . 0*AREAL( IELEM) * AREAL ( IELEM ) *ABS(DUDXA) 

EVALUATE RIGHT-HAND SIDE 

DO 50 1=1, NNODP 
K=L::ODS( IELEM. I) 

CARX I=DX ( I ) *DJA* CONSX*DTIME 
EQRHR ( K) =EQRHR f K ) -DDRDX*CARXI 
EQRHU ( K ) =EQRHU ( K ) -DRUDX *C ARX I 
EQRHE ( K ) =EQRHE ( K) -DREDX*CARXI 
50 CONTINUE 
100 CONTINUE 

RETURN 

END 


SUBROUTINE WRITER ( I ITER , RMSER , TSAVE) 

PARAMETER ( NELEM=200 , NP0IN=201 ) 

PARAMETER <NN0DP=2> 

COMMON / INIT/DENSY ( NPOIN ) , UVELY (NPOIN) , ENEGY (NPOIN) , PRESY (NPOIN) 
COMMON/COOR/XX ( NPOIN ) , LNODS ( NELEM , NNODP ) 

COMMON/ AAAA/A( NPOIN) ,DA(NPOIN) 

COMMON /TIME /CFLNB , DTI ME . ITMAX 
DIMENSION PB (NPOIN) , UB ( NPOIN) , RB ( NPOIN) 


WRITING PROCEDURES 


r* 


c 

c 


1010 


IF ( I ITER . EQ . 1 ) THEN 
READ (19,*) NUM 
I NUM= I TMAX / NUM 
END1F 

IF ( 7 ITER/ 200 *2 00 . EQ . I ITER ) 'WRITE ( 6 , 1000 ) I ITER , TSAVE , RMSER 

'WRITE (18, 1000 ) I ITER , TSAVE , RMSER 

NP1C 1=2 

NPIC2=NPOIN/2 

NPIC3=NPOIN-l 

NPIC2=23 

NPIC3=150 

NPIC4=200 

WRITE (18,1010) TSAVE, PRESY (NPIC1) . PRESY (NPIC2 ) , PRESY ( NPIC3 ) , 
& PRESY ( NP IC4 ) 

'WRITE (28, 1010) TSAVE, UVELY (NPIC1) . 

1 UVELY ( NPIC2 ) , UVELY(NPIC3 ) , UVELY ( NPIC4 ) 

FORMAT (5E12 .5) 

IF ( I ITER . EQ . 1 ) ICONT=INUM 
IF (I ITER. EQ. ICONT) THEN 
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ICONT=INUM 
DO 333 1=1 . NPOIN 
PB ( I ) =0 . 0 
UB ( I ) =0 . 0 
RB ( I ) =0 . 0 
TDIST=0 . 0 
333 CONTINUE 
END IF 

DO 444 1=1, NPOIN 
PB ( I ) =PB ( I ) +PRESY ( I ) *DTIME 
UB ( I ) =UB ( I ) +UVELY ( I ) *DTIME 
RB ( I ) =RB ( I ) +DENSY ( I ) *DTIME 
444 CONTINUE 

TD I ST=TD I ST+DT I ME 

IF( I ITER. EQ. ICONT) THEN 

DO 10 1=1, NPOIN 

PB ( I ) =PB ( I ) /TDIST 

UB ( I ) =UB( I ) /TDIST 

RB ( I ) =RB ( I ) ,/TDIST 

vmiTE( 16 , 1020) I , PB ( I ) , UB ( I) , RB ( I ) 

10 CONTINUE 

DO 15 1=1, NPOIN 
PB ( I ) =0 . 0 
UB ( I ) =0 . 0 
RB ( I ) =0 . 0 
15 CONTINUE 
TDIST=0 . 0 
I CONT= I I TER+ I NUM 
END IF 
C 

C WRITE AT EACH ICONT-TH ITERATION 

C 

C 

C ’WRITE IF SOLUTIONS ARE CONVERGED 
C 

IF (RMSER . GT . 1 . 0E-05 ) GO TO 20 
IF ( I ITER . EQ . 1 ) GO TO 20 
DO 30 1=1, NPOIN 

'WRITE (6 , 1020) I,XX(I) ,DENSY<1) ,UVELY(I> , ENEGY ( I ) ,PRESY(I) 
30 CONTINUE 

WRITE(13 , 1060) (XX(I) ,1=1, NPOIN) 

WRITE ( 13 , 1060) ( A ( I ) , 1=1, NPOIN) 

WRITE (13 , 1060) (DENSY ( I ) , 1=1 , NPOIN) 

WRITE (13 , 1060) ( UVELY ( I) , 1=1 , NPOIN) 

WRITE ( 13 , 1060) (ENEGY ( I ) ,1=1, NPOIN) 

WRITE (13 , 1060) ( PRESY ( I ) , 1=1 , NPOIN) 

STOP 

20 CONTINUE 
C 

C ’WRITE IF I ITER EQUALS TO ITMAX 
C 

IF ( I ITER. EQ. ITMAX) THEN 
DO 40 1=1, NPOIN 

C WRITE (6, 1020) I , XX ( I ) , DENSY ! I ) , UVELY ( I ) , ENEGY ( I ) , PRESY ( I ) 

40 CONTINUE 

'WRITE ( 13 . 1060) ( XX ( I ) , 1 = 1 , NPOIN) 

’WRITE ( 13 . 1060) <A< I ) , 1 = 1 .NPOIN) 

'WRITE (13. 1060) ( DENSY ( I ) , 1=1 .NPOIN) 

'WRITE (13. 1060) ( UVELY ( I ) , 1 = 1. NPOIN) 

WRITE ( 13 , 1060) (ENEGY ( I ) , 1=1 , NPOIN) 


WRITE ( 13 , 1060 ) (PRESY(I) ,1=1, NPOIN) 
END IE 
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C 

RETURN 

1000 FORMAT (5X, 15 , 2X, (2 (E10 . 5 , IX) ) ) 

1020 FORMAT (5X,I5,F10.5,4E15.5) 

1060 FORMAT ( 5 ( 200 ( 4E15 . 5 , / ) ) ) 

END 

C 

C 

SUBROUTINE ITERAT 

PARAMETER ( NELEM=200 , NP0IN=201 ) 

PARAMETER ( NN0DP=2 , NEQNS=3 ) 

PARAMETER (NTIME=2000) 

C PARAMETER (NTHNN=20 , NTHEE=NTHNN-1 ) 

PARAMETER ( NTHNN=201 , NTHEE=NTHNN-1 ) 

COMMON/MASS/GMASS (NPOIN) 

COMMON/AAAA/A( NPOIN) , DA (NPOIN) 

COMMON /DOMA/ DO , UO , EO , PO 
COMMON /BCBC/D1 , U1 , El , PI 
COMMON/ AREA/ AREAL ( NELEM) 

COMMON / INIT/DENSY ( NPOIN ) , UVELY(NPOIN) , ENEGY (NPOIN) , PRESY(NPOIN) 
COMMON/COOR/XX( NPOIN) . LNODS ( NELEM, NNODP) 

COMMON/ PRTY/CAPAV , CAPAP , CGAM , CONDT , VI SCY 
COMMON/TIME/CFLNB , DTIME , ITMAX 

COMMON /EQNS/EQRHR< NPOIN) , EQRHU ( NPOIN ) , EQRHE (NPOIN ) 

DIMENSION DELTR ( NPOIN ) , DELTU ( NPOIN) , DELTE ( NPOIN ) , 

DENST ( NPOIN ) , UVELT(NPOIN ) , ENEGT (NPOIN) 

DIMENSION EQRHS(NPOIN) 

DIMENSION DDS ( NTIME , NTHNN ) , PPS ( NTIME , NTHNN) , UUS ( NTIME , NTHNN ) 
DIMENSION NNSS(NTHEE, NNODP) ,XXS(NTHNN) , AS ( NTHNN) 

DIMENSION PSTA( NTHNN) , USTA ( NTHNN) , DSTEP ( NTIME ) 

DIMENSION PBAR( NTHNN) ,UBAR( NTHNN) ,RBAR( NTHNN) 

DIMENSION RSTA( NTHNN) 

C 

OPEN ( 16 , FILE=' Stl6.dat' ) 

NELS=NTHEE 

NXS=NTHNN 

NTS=NTIME 

DO 441 1=1 , NELS 

DO 441 J=1 , 2 

NNSS ( I , J ) =LNODS ( I , J ) 

441 CONTINUE 

DO 442 1=1, NXS 
XXS(I)=XX(I) 

AS ( I ) =A ( I ) 

PSTA( I ) =0 . 

USTA ( I ) =0 . 

RSTA ( I ) =0 . 

442 CONTINUE 
RR1=DENSY ( 1 ) 

UU1=UVELY ( 1 ) 

EE1=ENEGY ( 1 ) 

PP1=PRESY ( 1 ) 

TT1=PP1/RR1/ ( CGAM-1 . ) /CAPAV 
ASOUND=SQRT (CGAM*PP1/RR1) 

AMACH=UU1/ ASOUND 
THLEN=XX ( 20 ) -XX ( 1 ) 

C THLEN=XX ( NPOIN ) -XX ( 1 ) 

FREQ=3 . 14* ASOUND/' THLEN 


PRINT* , ASOUND , AMACH , FREQ 

C 

C SET UP ITERATION COUNTER AND LOOP ADRESS 
C 

READ (19,*) CPERT 
I ITER=0 
IC0UN=0 
C 

DDTM=0 . 

TSAVE=0 . 0 
10 CONTINUE 

I ITER=I ITER+1 

C 

C SET UP TIME STEP (VARIABLE TIME STEP) 

C 

C IF (I ITER. EQ. 1) TSAVE=DTIME 

IF (I ITER. GE. 1) CALL SDTIME 
C DTIME=1 . OE-5 

IF (I ITER. GE. 1) T$AVE=TSAVE+DTIME 

C 

RCOS=SIN ( FREQ*TSAVE ) 

C RC0S=1 . 0 

PERT=CPERT* PP1 
NTRIC=1000 

IF ( I ITER . GE . NTRIG) PERT=CPERT* PP1 

TPRE=PP1+PERT*RC0S 

CIGAM=1 . /CGAM 

PRRR=TPRE/PP1 

CC PRINT* , I , RRAD, RCOS .TPRE 
TTEM=TT 1/PRRR* *l’I GAM 
TUUU=UU1 

TSQR=0 . 5*TUUU**2 
TENG=CAPAV*TTEM+TSQR 
TRHO=TPRE/ ( ( CGAM-1 . ) * (TENG-TSQR) ) 
TENT=TENG+TPRE/TRHO 
DENSY ( 1 ) =TRHO 
UVELY ( 1 ) =TUUU 
ENEGY ( 1 ) =TENG 
PRESY ( 1 ) =TPRE 
C 

C INITIALIZATIONS 
C 

C GOTO (21.22,23), IEQNS 

21 DO 25 1=1 , NPOIN 

25 EQRHR ( I ) =0 . 0 

C GO TO 20 

22 DO 26 1=1, NPOIN 

26 EQRHU ( I ) =0 . 0 

C GO TO 20 

^ 23 DO 27 1=1 , NPOIN 

27 EQRHE ( I ) =0 . 0 
20 CONTINUE 

C 

C ASSEMBLE CONTRIBUTIONS OF EACH ELEMENT TO THE RIGHT-HAND 
C SIDE VECTORS 
C 
C 

C DOMAIN CONTRIBUTIONS 
C 

DO 30 IELEM=1 , NELEM 


o o o ooo ooo 


75 


CALL MATRIX ( I ITER . 1 , IELEM) 

30 CONTINUE 

SURFACE CONTRIBUTIONS 

CALL BDFLUX(l) 

MAIN MATRIX SOLVER USING ITERATIVE SCHEME 
DO 90 IEQNS=1 , NEQNS 

IF ( IEQNS . EQ . 1 ) CALL SOLVER ( I EQNS , DELTR) 

IF ( IEQNS . EQ . 2 ) CALL SOLVER ( IEQNS , DELTU) 

IF ( IEQNS . EQ. 3 ) CALL SOLVER ( IEQNS, DELTE) 

90 CONTINUE 

UPDATE SOLUTIONS 

DO 100 1=1 , NPOIN 
DENST ( I ) =DENSY ( I ) *A ( I ) +DELTR ( I ) 

UVELT ( I ) =DENSY ( I ) *UVELY ( I ) *A ( I ) +DELTU ( I ) 

ENEGT ( I ) =DENSY ( I ) * ENEGY ( I ) *A ( I ) +DELTE ( I ) 

100 CONTINUE 

DO 110 1=1, NPOIN 
DENSY ( I ) =DENST ( I ) / A ( I ) 

UVELY ( I ) = UVELT ( I ) /DENST ( I ) 

ENEGY ( I ) =ENEGT ( I ) /DENST ( I ) 

PRESY ( I ) = (CGAM-1 . 0 ) *DENSY ( I ) * ( ENEGY ( I ) -0 . 5 *UVELY ( I ) * *2 ) 
C PRINT* , I , DENSY ( I ) , UVELY ( I ) , ENEGY ( I ) , PRESY ( I ) 

110 CONTINUE 

C 

C CHECK THE CONVERGENCE 

C 

SUMUP=0 . 0 
SUMDN=0 . 0 
DO 115 1=1, NPOIN 

SUMUP=SUMUP+DELTR ( I ) * * 2 +DELTU ( I ) * * 2 +DELTE i I ) * * 2 
SUMDN=SUMDN+DENST ( I ) * *2+UVELT ( I ) * *2+ENEGT ( I ) **2 
115 CONTINUE 

RMSER=SQRT ( SUMUP / SUMDN ) 

C 

C APPLY LAPIDUS' ARTIFICIAL VISCOSITY 
C 

C DO 190 IEQNS=1, NEQNS 

C GOTO (121,122,123), IEQNS 

121 DO 125 1=1, NPOIN 

125 EQRHR ( I ) =0 . 0 

C GO TO 120 

122 DO 126 1=1. NPOIN 

126 EQRHU ( I ) =0 . 0 

C GO TO 120 

123 DO 127 1=1, NPOIN 

127 EQRHE ( I ) =0 . 0 
120 CONTINUE 

C 

CALL LAPDUS ( 1 ) 

C GO TO (140.150,160). IEQNS 

C 

140 DO 145 1=1, NPOIN 

145 DENST ( I ) =DENSY ( I ) * A ( I ) +EQRHR ( I ) /GMASS ( I ) 

C GO TO 180 
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150 DO 155 1=1 . NPOIN 

155 UVELT ( I ) =DENSY ( I ) * UVELY ( I ) *A ( I ) +EQRHU ( I ) /GMASS ( I ) 

GO TO 180 

160 DO 165 1=1, NPOIN 

165 ENEGT ( I ) =DENSY ( I ) * ENEGY ( I ) * A ( I ) +EQRHE ( I ) /GMASS ( I ) 

180 CONTINUE 
190 CONTINUE 

COMPUTE FINAL SOLUTIONS AT EACH TIME STEP 

DO 200 1=1, NPOIN 
DENSY ( I ) =DENST ( I ) /A ( I ) 

UVELY ( I ) =UVELT ( I ) /DENST ( I ) 

ENEGY ( I ) =ENEGT ( I ) /DENST ( I ) 

PRESY ( I ) = ( CGAM-1 , 0 ) *DENSY ( I ) * ( ENEGY ( I ) -0 . 5 *UVELY ( I ) * * 2 
200 CONTINUE 

SUB SON i.O INLET BOUNDARY CONDITIONS 

CASE A 

Dl=DO 

U1-U0 

El=EO 

P1=°0 

CASE C 

HO= (CGAM/ (CGAM-1 . 0 ) ) *PO/D0+0 . 5*U0*U0 
CGAM1=CGAM-1 . 0 
R3 ^0/D0**CGAM 

D1 - ( (CGAM1/CGAM) * (HO-0 . 5*UVELY (1)**2)/R3)**(1. /CGAM1 ) 
Ul= JVELY(l) 

P1=.13*D1**CGAM 

E1--’I/ ( (CGAM-1 .0) *Dl)+0. 5*U1*U1 

DENSY ( 1 ) =D1 

UVELY(1)=U1 

PRESY ( 1 ) =P1 

ENEGY ( 1 ) =E1 

SUBSONIC OUTLET BOUNDARY CONDITIONS 

PRESY ( NPOIN )=0. 704 
I ' 3 RESY (NPOIN) =0 , 61845265 

CALL WR [TER TO OUTPUT ITERATION RESULTS 

C 

CAL^ ’WRITER ( I ITER , RMSER , TSAVE ) 

DD n.I=DDTM+DTIME 
NONE=ITMAX/2000 
IA=I ITER/NONE 
IB=NONE* IA 

IF ( I ITER . EQ . IB ) THEN 
I COUN = I COUN + 1 
DSTEP ( ICOUN ) =DDTM 
DDTM=0 . 

DO 957 1=1, NXS 

PPS( ICOUN, I ) =PRESY ( I ) 

UU3 ( ICOUN , I ) = UVELY ( I ) 

DDS ( ICOUN . I ) =DENSY ( I ) 

957 CONTINUE 
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443 

C 


871 

C 


C 


c 


c 


c 


c 


c 

c 


c 


c 

c 

c 


ENDIF 

DO 443 1=1, NXS 

PSTA ( I ) = PSTA ( I ) +DTIME*PRESY ( I ) 

USTA ( I ) =USTA ( I ) +DTIME*UVELY ( I ) 

RSTA ( I ) =RSTA ( I ) +DTIME*DENSY ( I ) 

CONTINUE 

IF ( IITER.LT. ITMAX) GO TO 10 
DO 871 1=1, NXS 
PSTA ( I ) =PSTA ( I ) /TSAVE 
USTA ( I ) =USTA ( I ) /TSAVE 
RSTA ( I ) =RSTA ( I ) /TSAVE 
CONTINUE 
REWIND (16) 

CALL STAB ( NELS , NXS t NTS , XXS , AS , PPS , UUS , DDS , NNSS , DSTEP , 
PSTA , USTA , RSTA , PBAR , UBAR , RBAR ) 

RETURN 

END 

SUBROUTINE STAB ( NELE , NX , NT , X , A , P , U , R , NEL t DSTEP , 

PSTA , USTA , RSTA , PBAR , UBAR , RBAR ) 


PARAMETER ( NINT=2 , L=2 , NF=12 ) 

COMMON/ PARAMT/ SOUND , GAMMA , ADM I , ADMO 

COMMON/ PVALUE/XLENG , TLENG , I PERT , PI , EPSI , VISCO 

DIMENSION NEL (NELE, L) 

DIMENSION X ( NX ) , A ( NX ) , P ( NT , NX ) , U ( NT , NX ) , R ( NT , NX ) 
DIMENSION XI (NINT) ,WI (NINT) , PHI ( L , NINT) 

DIMENSION DSTEP(NT) ,PSTA(NXj RSTA ( NX ) , USTA ( NX ) 

DIMENSION PBAR ( NX ), UBAR (NX) RBAR (NX) 

IPERT=0 

PI=3 . 141592654 
GAMMA=1 . 2 
ADMI =0 . 0 
ADMO=0 . 0 
VISC0=0. 0 
EPSI=0. 2 

XLENG=X ( NX ) -X ( 1 ) 

TLENG=0 . 05 

CALL GAUSS (NINT, XI ,WI) 

CALL SHAPE ( NINT , XI , PHI ) 

CALL ABC ( NELE , L , NF , NX , NINT , NT , NEL , X , A , R , P , U , 

WI , PH I , AA , BB , CC , DSTEP , PSTA , USTA , RSTA , PBAR , UBAR , RBAR) 

RETURN 

END 


SUBROUTINE ABC (NELE , L , NF , NX , NINT , NT , NEL , X , A , R , P , U , WI , PHI , AA , BB , CC , 
DSTEP . PSTA , USTA , RSTA , PBAR , UBAR , RBAR ) 

COMMON /PARAMT/ SOUND , GAMMA , ADMI , ADMO 
COMMON/PVALUE/XLENG, TLENG, IPERT, PI , EPSI , VISCO 


C 



DIMENSION NEL (NELE , L) ,WI(NINT) ,PHI(L,NINT) 
DIMENSION X ( NX ) ,A(NX) ,P(NT,NX) ,U(NT,NX) ,R(NT,NX) 
DIMENSION E ( 3 ) ,G(3) ,ETA(3) ,EE<3) ,GG(3) ,ET(3) 

C 

DIMENSION PBAR(NX) ,UBAR(NX) , RBAR(NX) , DSTEP (NT) 
DIMENSION PSTA(NX) ,RSTA(NX) ,USTA(NX) 

DIMENSION ICAL( 10000) , IDAL( 10000) 

C 

NXS=NX 

NTS=NT 

C 

READ( 19 , * ) INUM 
NUM=2000/INUM 
DO 941 1=1, NT 
ICAL ( I ) =0 
IA=I/NUM 
IB=NUM* IA 
IF(I.EQ.IB) THEN 
ICAL ( I ) =1 
END IF 

941 CONTINUE 
ICAL ( 1 ) =1 
ICAL (NT) =0 

C 

AAA=0.0 
BBB=0 . 0 
CCC=0.0 
ITER=0 

101 ITER=ITER+1 

PRINT* ITRT= , ITER 
DO 102 1=1,3 
EE ( I ) =0 . 0 
GG ( I ) =0 . 0 
ET ( I ) =0 . 0 
102 CONTINUE 
C 
C 

C NT=100 

C STEP=TLENG/NT 

C DO 200 1=1 , NT 

C T= ( 1-1 ) *STEP 

C PRINT*, (DSTEP(I) , 1=1, NT) 

S0UND=0 . 0 
RHHH=0 . 0 
EZZZ=0 . 0 
C 

T=0 . 0 
TZERO=T 
NTRIG=600 • 

^ DO 100 I ITER=1 , NT 

IF ( ICAL ( I ITER) . EQ . 1 ) THEN 
DO 660 1=1, NX 

READ (16,*) NUM , PSTA ( I ) , USTA ( I ) , RSTA ( I ) 

660 CONTINUE 
END IF 
XSND=0 . 0 
,KRRR=0 . 0 

xyuu=o.o 

DO 443 1=1 . NXS 
PBAR ( I ) =P ( I ITER , I ) 
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UBAR ( I ) =U ( I ITER , I) 

RBAR ( I ) =R ( I ITER , I) 

XRRR=XRRR+PBAR ( I ) 

XUUU=XUUU+UBAR ( I ) 

XSND=XSND+SQRT ( GAMMA* PB AR ( I ) / RBAR ( I ) ) 
443 CONTINUE 

SOUND=SOUND+XSND/FLOAT ( NXS*NTS) 
RHHH=RHHH+XRRR/ FLOAT ( NXS * NTS ) 
SOUND=XSND/ FLOAT ( NXS ) 

RHHH=XRRR/ FLOAT ( NXS ) 

XPPP=XRRR/ FLOAT ( NXS ) 

PRINT* , SOUND , RHHH , XPPP 

XUUU=XUUU/ FLOAT ( NXS ) 

STEP=DSTEP ( I ITER ) 

T=T+STEP 


CALL DOMAIN (NELE , L , NINT , NEL , NX , 'JVI . PHI ,T,X, A, 

- PBAR , UBAR , RBAR , E , G , PSTA , USTA , RSTA) 

C 

DO 250 J=1 , 3 

EE ( J ) =EE ( J ) +STEP*E ( J ) 

GG ( J ) =GG ( J ) +STEP*G ( J ) 

250 CONTINUE 
C 

CALL BOUND ( NX, T,X, A, PBAR, UBAR. RBAR, ETA, PSTA, USTA, RSTA) 

C 

DO 350 J=1 . 3 

ET ( J ) =ET ( J ) +STiP*ETA ( J ) 

350 CONTINUE 
C 

IF ( ICAL ( I ITER) . EQ . 1 ) THEN 
ONE=ET ( 1 ) +GG ( 1 ) 

TWO=ET ( 2 ) +GG ( 2 ) 

THR=ET ( 3 ) +GG ( 3 ) 

C 

CONE=0 . 5 / EE ( 1 ) 

CTW0=EE(2)/EE( I) 

CTHR=EE ( 3 ) /EE ( I ) 

AA=ONE*CONE 

BB= (TWO-1 . 5*0NE*CTW0) *CONE 

CC= ( THR-1 , 5 *TWO*CTWO+ (2.2 5 *CTWO* * 2-2 . 0 *CTHR ) *ONE ) *CONE 
WRITE ( 24 , 1031) T,AA,BB,CC 
C CALL EULER ( AA , BB . CC , TZERO , TEND) 

DO 107 I J=1 . 3 
EE ( I J ) =0 . 0 
GG ( I J ) =0 . 0 
ET ( I J ) =0 0 
107 CONTINUE 
END IF 

100 CONTINUE 
1031 FORMAT ( 2X, 4E14 5) 

C 

C PRINT* . I ITER . EE ( 1 ) , EE ( 2 ) , EE ( 3 ) 

PRINT* . EE ( 1 ) , EE ( 2 ) , EE ( 3 ) 

PRINT* .CONE, CTWO.CTHR 
PRINT* , AA, BB , CC 
TEND=T 
C 
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DDA=ABS ( AA-AAA ) 

DDB=ABS ( BB-BBB ) 

DDC=ABS (CC-CCC ) 

RMA=AA* *2+BB**2+CC**2 

RMB=DDA * * 2 +DDB * * 2 +DDC * * 2 
RMC=SQRT ( RMB/RMA) 

AAA=AA 
BBB=BB 
CCC=CC 
TZERO=TEND 

PRINT* , ITER , DDA , DDB , DDC , RMC 
IF ( RMC . GT . 1 . OE-4 . AND . ITER . LT . 1 ) GO TO 101 
PRINT*, ITER, RMC 
RETURN 
END 
C 
C 
C 
C 

SUBROUTINE EULER ( AA , BB , CC , TZERO , TEND ) 

C 

COMMON/ PARAMT/ SOUND , GAMMA , ADMI . ADMO 

COMMON/ PVALUE/XLENG , TLENG , I PERT, PI , EPS I .VI SCO 


NT=100 
EPSI=1 . 0 

STEP= (TEND-TZERO ) / (NT-1> 

DO 100 1=1, NT 
T=TZERO+ ( I - 1 ) * STEP 
EZERO=EPSI 
EEND=0 . 0 
1000 CONTINUE 

C0NE=STEP/2 . 0* ( AA+2 . 0*BB*EEND+3 . 0*EEND* *2 ) 
CTW0=EEND-STEP/2 . 0* ( AA*EEND+BB*EEND+CC*EEND* *3 ) 
CTHR=-EZER0-STEP/2 . 0* ( AA*EZERO+BB*EZERO* *2+CC*EZER0* *3 ) 
DELTE=- ( CTWO+CTHR ) / ( 1 . O-CONE ) 

EEND=EEND+DELTE 

IF(ABSfDELTE) . LT.l.OE— 5) GOTO 1000 
WRITE (26, 11) T.EEND 
11 FORMAT ( 2E15 . 3 ) 

100 CONTINUE 


C 


C 


C 


C 


RETURN 

END 


SUBROUTINE DOMAIN ( NELE , L , NINT , NEL , NX , WI , PHI , T , X , A , P , U . R , E , G , 
PSTA , USTA , RSTA ) 

COMMON/PARAMT/SOUND , GAMMA , ADMI , ADMO 

COMMON/ PVALUE/XLENG, TLENG, IPERT, PI , EPSI , VISCO 

DIMENSION NEL (NELE, L) , X (NX > , A< NX ) , P < NX ) , U< NX > , R < NX > 

DIMENSION WI ( NINT) ,PHI(L,NINT) , DPDS ( 2 ) , EE ( 5 ) 

DIMENSION QX ( 201 ) ,QA(201) ,QP(201) ,QU(201) ,QR(201) 

DIMENSION QPR ( 201 ) ,QUR(201) ,QRR(201) 

DIMENSION E ( 3 ) , G ( 3 ) , ENT ( 5 ) , DE ( 201 , 5 ) 

DIMENSION PSTA (NX) , RSTA (NX) , USTA (NX) ,S(4) 

DIMENSION RV ( 5 ) ,DEX(5) , DPVX ( 5 ) , PR ( 5 ) , PRV ( 5 ) ,DRX(5) 


DIMENSION RW ( 5 ) , DVX ( 5 ) 


DPDS ( 1 ) =-0 . 5 
DPDS ( 2 ) =0 . 5 
DO 50 1=1,3 
E(I) =0.0 
G ( I ) =0 . 0 
50 CONTINUE 

DO 100 1=1, NX 
QX ( I ) =X ( I ) 

QA ( I ) =A ( I ) 

QP ( I ) =PSTA ( I ) 

QPR ( I ) =P ( I ) -PSTA ( I ) 

QU ( I ) =USTA ( I ) 

QUR ( I ) =U ( I ) -USTA ( I ) 

QR(I)=RSTA( I) 

QRR ( I ) =R ( I ) -RSTA ( I ) 

PP=QP( I ) 

P1’R=QPR( I ) 

UU=QU(I) 

UPR=QUR ( I ) 

RR=QR ( I ) 

RPR=QRR ( I ) 

C^LL ENTHAL (PP , PPR , UU , UPR , RR , RPR , EE ) 
DO 150 J=1 , 5 
DE ( I , J ) =EE ( J ) 

150 CONTINUE 
100 CONTINUE 

DO 200 1=1 , NELE 
DO 250 J=1 , NINT 

x.:=o.o 

i \ ji = 0 . 0 

PP=0. 0 

p:>r=o. o 

uu=o . 0 

LVR=0 . 0 
R)t=0 . 0 
PPR=0 . 0 

DO 280 JONE=l , 5 
ENT ( JONE) =0 . 0 
280 CONTINUE 

DO 300 K=1,L 
.'I ,T M=NEL ( I , K) 

CON=PHI (K., J) 

XX=XX+QX ( NUM ) *C0N 
A\=AA+QA(NUM) »C0N 
UU=UU+QU ( NUM ) *CON 
UPR=UPR+QUR ( NUM ) *CON 
PP=PP+QP ( NUM ) *C0N 
PPR=PPR+QPR ( NUM ) *C0N 
RR=RR+QR ( NUM) *CON 
RPR=RPR+QRR( NUM) *C0N 
DO 350 K0NE=1,5 

ENT ( KONE ) =ENT ( KONE ) +DE ( NUM , KONE ) *CON 
350 CONTINUE 
300 CONTINUE 
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C JCB = -0 . 5 •* X ( NEL (1,1) ) +0 . 5*X (NEL (1.2) ) 

C 

CONE=CJCB*WI ( J) * AA 
C 

E ( 1 ) =E ( 1 ) + ( RR*ENT ( 3 ) +RPR*ENT ( 2 ) ) *CONE 
E ( 2 ) =E ( 2 ) + ( RR*ENT ( 4 ) +RPR*ENT ( 3 ) ) *CONE 
E ( 3 ) =E ( 3 ) + ( RR*ENT ( 5 ) +RPR*ENT ( 4 ) ) *C0NE 
C 

DO 400 KONE-1,5 
DEX ( KONE ) =0 . 0 
DO 450 KTWO=l,L 
NUM=NEL ( I , KTWO ) 

DEX ( KONE ) =DEX ( KONE ) +DE ( NUM , KONE ) * DPDS ( KTWO ) 

450 CONTINUE 

DEX ( KONE ) =DEX ( KONE ) /C JCB 
400 CONTINUE 
C 

DO 480 K=1 , 3 
DPVX (K) =0 . 0 
480 CONTINUE 
C 

DO 550 KTW0=1.L 
NUM=NEL ( I . KTWO ) 

DPVX ( 1 ) =DPVX ( 1 ) +QP (NUM) *QU(NUM) *DPDS (KTWO) 

DPVX ( 2 ) =DPVX ( 2 ) + ( QP ( NUM ) *QUR ( NUM) +QPR < NUM ) *QU(NUM) ) *DPDS(KTWO) 
DPVX ( 3 ) =DPVX ( 3 ) +QPR ( NUM ) *QUR ( NUM ) * DPDS ( KTWO ) 

550 CONTINUE 
C 

DO 600 K0NE=1,3 
DPVX ( KONE ) =DPVX ( KONE ) /C JCB 
600 CONTINUE 
C 

DO 650 K0NE=1 , 2 
DRX( KONE) =0.0 
650 CONTINUE 

DO 700 K0NE=1,2 
NUM=NEL ( I , KONE ) 

DRX ( 1 ) =DRX ( 1 ) +QR ( NUM ) *DPDS ( KONE ) 

DRX ( 2 ) =DRX ( 2 ) +QRR ( NUM ) * DPDS ( KONE ) 

700 CONTINUE 

DO 750 K0NE=1 , 2 
DRX ( KONE ) =DRX ( KONE ) /C JCB 
750 CONTINUE 
C 

DO 810 K0NE=1 , 2 
DVX ( KONE ) =0 . 0 
810 CONTINUE 
C 

DO 820 K0NE=1, 2 
^ NUM=NEL ( I . KONE ) 

DVX ( 1 ) =DVX ( 1 ) +QU ( NUM ) *DPDS (KONE) 

DVX ( 2 ) =DVX ( 2 ) +QUR ( NUM ) *DPDS ( KONE ) 

820 CONTINUE 
C 

DO 830 K0NE=1 , 2 
DVX ( KONE ) =DVX ( KONE ) /C JCB 
830 ,'GONT I NUE 
C 

RV ( 1 ) =RR*UU 

RV ( 2 ) =RR*UPR+RPR*UU 


EV ( 3 ) =RPR*UPR 

C 

CA=1.0/( GAMMA-1.0) 

CB=-GAMMA/ ( GAMMA-1 . 0 ) 

PC=PPR/PP 
RC=RPR/RR 
S ( 1 ) =CA* PC+CB*RC 

S(2)=-l. 0/2.0* (CA*PC**2+CB*RC**2 ) 

S(3) =10/6.0* (CA*PC**3+CB*RC* *3 ) 

S ( 4 ) =-l .0/48.0* (CA*PC**4+CB*RC**4 ) 

C 

C0N=1 . 0/RR**4 
PR< 1 ) =PP*RR* *3 *CON 

PR < 2 ) = ( -PP*RR* * 2 * RPR+RR* * 3 * PPR ) *C0N 
PRO) = (PP*RR*RPR**2-RR**2*PPR*RPR) *C0N 
PR(4) = ( -PP*RPR* *3+RR*PPR*RPR* *2 ) *C0N 
PR(5) =-PPR*RPR**3*C0N 

C 

PRV ( 1 ) =PR ( 1 ) *UU 
PRV ( 2 ) =PR ( 1 ) *UPR+PR ( 2 ) *UU 
PRV ( 3 ) =PR ( 2 ) *UPR+PR ( 3 ) *UU 
PRV ( 4 ) =PR ( 3 ) *UPR+PR ( 4 ) *UU 
PRV ( 5 ) =PR ( 4 ) *UPR+PR ( 5 ) *UU 
C 

RW ( 1 ) =RR*UU* *2 

RW ( 2 ) =2 . 0*RR*UU*UPR+UU* * 2 *RPR 
RW ( 3 ) =RR*UPR**2+2 . 0*UU*UPR*RPR 
RW ( 4 ) =RPR*UPR* * 2 
C 

G ( 1 ) = G ( 1 ) + ( RV ( 1 ) *DEX ( 3 ) +RV ( 2 ) *DEX ( 2 ) +RV ( 3 ) *DEX ( 1 , ) *C0NE 
G(2 ) =G( 2 ) + ( RV ( 1 ) *DEX ( 4 ) +RV ( 2 ) *DEX ( 3 ) +RV ( 3 ) *DEX (2 > ) *C0NE 
G ( 3 ) =G ( 3 ) + ( RV ( 1 ) *DEX ( 5 ) +RV ( 2 ) *DEX ( 4 ) +RV ( 3 ) *DEX ( 3 ) ) *C0NE 
C 

G ( 1 ) =G < 1 ) + ( DPVX ( 1 ) * S ( 2 ) +DPVX ( 2 ) * S ( 1 ) ) *C0NE 
G ( 2 ) =G ( 2 ) + { DPVX ( 1 ) * S ( 3 ) +DPVX ( 2 ) * S ( 2 ) +DPVX ( 3 ) *S ( 1 ) ) *C0NE 
G ( 3 ) =G ( 3 ) + ( DPVX ( 1 ) * S ( 4 ) +DPVX ( 2 ) * S ( 3 ) +DPVX ( 3 ) * S ( 2 ) ) *C0NE 
C 

G( 1 ) =G( 1 ) -(DRX ( 1 ) *PRV ( 3 ) +DRX ( 2 ) *PRV ( 2 ) ) *C0NE 
G ( 2 ) =G ( 2 ) - (DRX ( 1 ) *PRV ( 4 ) +DRX ( 2 ) *PRV ( 3 ) ) *C0NE 
G ( 3 ) =G ( 3 ) - ( DRX ( 1 ) * PRV ( 5 ) +DRX ( 2 ) * PRV ( 4 ) ) *C0NE 

C 

G ( 1 ) =G ( 1 ) - ( DVX ( 1 ) *RW ( 3 ) +DVX ( 2 ) * RW ( 2 ) ) *C0NE 
G ( 2 ) =G ( 2 ) - ( DVX ( 1 ) * RW ( 4 ) +DVX ( 2 ) * RW ( 3 ) ) *C0NE 
G( 3 ) =G ( 3 ) -DVX ( 2 ) *RW ( 4 ) *C0NE 
C 

250 CONTINUE 
200 CONTINUE 

C 

RETURN 

END 

C 

C 

C 

SUBROUTINE BOUND (NX ,T,X,A,P,U,R, ETA , PSTA, USTA , RSTA ) 

C 

COMMON/ PARAMT/ SOUND , GAMMA . ADM I , ADMO 
COMMON /PVALUE/XLENG , TLENG , IPERT, PI , EPSI ,VISCO 
C 

DIMENSION X ( NX ) , A (NX) , P (NX) ,U(NX) , R(NX) 

DIMENSION EE ( 5 ) , ETA( 3 ) 
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DIMENSION PSTA(NX) ,RSTA(NX) ,USTA(NX) ,S(4) 
DIMENSION RV ( 5 ) , PV ( 5 ) 

C 

DO 1201 1=1,3 
ETA ( I ) =0 . 0 
1201 CONTINUE 
E0NE=0 . 0 
ETWO=0 . 0 
ETHR=0 . 0 
C 

XX=X(1) 

PP=PSTA ( 1 ) 

UU=USTA ( 1 ) 

AA=A ( 1 ) 

RR=RSTA ( 1 ) 

C 

PPR=P ( 1 ) -PSTA ( 1 ) 

UPR=U ( 1 ) -USTA ( 1 ) 

RPR=R ( 1 ) -FfoTA ( 1 ) 

C 

CALL ENTHAL ( PP . PPR , UU , UPR , RR , RPR , EE ) 

C 

RV ( 1 ) =RR*UU 
RV ( 2 ) =RR*UPR+RPR*UU 
RV ( 3 ) =RPR*UPR 
C 

CONE=RV ( 1 ) *EE ( 3 ) +RV ( 2 ) *EE ( 2 ) +RV ( 3 ) *EE ( 1 ) 
CTWO=RV ( 1 ) *EE ( 4 ) +RV ( 2 ) *EE(3)+RV(3) *EE(2) 
CTHR=RV ( 1 ) "EE ( 5 ) +RV ( 2 ) *EE(4)+RV(3) *EE(3) 

C 

EONE=EONE-CONE 
ETWO=ETWO-l TWO 
ETHR=ETHR-CTHR 
C 

PV(1)=PP*UU 
PV ( 2 ) =PP*UPR+PPR*UU 
PV ( 3 ) =PPR*UPR 
C 

CA=1 .0/ (GAJ.ftlA-1 .0) 

CB=-GAMMA/ (GAMMA-1 .0) 

PC=PPR/PP 
RC=RPR/RR 
S ( 1 ) =CA*PC CB*RC 

S(2)=-l. 0/2.0* (CA*PC**2+CB*RC**2 ) 

S ( 3 ) =1 . 0/6 0* (CA*PC**3+CB*RC**3 ) 
S(4)=-1.0/',8.0*(CA*PC**4+CB*RC**4) 

C 

EONE=EONE- i PV ( 1 ) *S(2)+PV(2) *S(1) )-PPR*UPR 
ETWO=ETWO- ' PV ( 1 ) * S ( 3 ) +PV ( 2 ) *S(2 ) +PV (3) *S( 1) ) 
ETHR=ETHR- ( PV ( 1 ) *S (4 ) +PV (2 ) *S (3 ) +PV (3 ) *S ( 2 ) ) 
C 

ETA ( 1 ) =EONE* AA 
ETA ( 2 ) =ETW0* AA 
ETA ( 3 ) =ETHR*AA 

C 

E0NE=0 . 0 
ETW0=0 . 0 
ETHR=0 . 0 
C 

XX=X ( NX ) 
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PP=PSTA(NX) 

UU=USTA ( NX ) 

AA=A (NX ) 

RR=RSTA(NX) 

PPR=P ( NX ) -PSTA ( NX ) 

UPR=U ( NX) -USTA(NX) 

RPR=R ( NX ) -RSTA ( NX ) 

C 

CALL ENTHAL ( PP , PPR , UU , UPR , RR , RPR , EE ) 

C 

RV ( 1 ) =RR*UU 
RV(2) =RR*UPR+RPR*UU 
RV ( 3 ) =RPR*UPR 
C 

CONE=RV ( 1 ) *EE ( 3 ) +RV ( 2 ) *EE ( 2 ) +RV ( 3 ) *EE ( 1 ) 

CTWO=RV ( 1 ) *EE ( 4 ) +RV ( 2 ) *EE ( 3 ) +RV ( 3 ) *EE ( 2 ) 

CTHR=RV(1) *EE (5 ) +RV(2) *EE ( 4 ) +RV ( 3 ) *EE (3) 

C 

EONE=EONE-CONE 

ETWO=ETWO-CTWO 

ETHR=ETHR-CTHR 

C 

PV(1 ) =PP*UU 

PV ( 2 ) =PP*UPR+PPR*UU 

PV( 3 ) =PPR*UPR 

C 

CA=1 . 0/ (GAMMA-1 .0) 

CB=-GAMMA/( GAMMA-1. 0) 

PC=PPR/PP 

RC=RPR/RR 

S ( 1 ) =CA*PC+CB*RC 

S ( 2 ) =-0 . 5* (CA*PC* *2+CB*RC* *2 ) 

S(3)=l . 0/6 . 0* (CA*PC**3+CB*RC**3 ) 

S(4) =-1.0/48.0* (CA*PC* *4+CB*RC* *4 ) 

C 

E0NE=E0NE- ( PV ( 1 ) * S ( 2 ) +PV ( 2 ) * S ( 1 ) ) -PPR*UPR 
ETW0=ETW0- < PV ( 1 ) * S ( 3 > +PV ( 2 ) * S ( 2 ) +PV ( 3 ) *S ( 1 ) ) 
ETHR=ETHR- (PV(1)*S(4)+PV(2)*S(3)+PV(3)*S(2)) 

C 

ETA ( 1 ) =-ETA ( 1 ) +E0NE* AA 
ETA ( 2 ) =-ETA ( 2 ) +ETW0* AA 
ETA ( 3 ) =-ETA ( 3 ) +ETHR* AA 
C 

RETURN 

END 

SUBROUTINE ENTHAL ( PP . PPR , UU . UPR , RR . RPR . EE ) 

C 

COMMON/ PARAMT/ SOUND , GAMMA , ADMI , ADMO 
COMMON/PVALUE/XLENG , TLENG , I PERT, PI , EPSI , VISCO 

C 

DIMENSION EE ( 5 ) 

C 

CON= (GAMMA-1 .0) * RR* *4 

EE ( 1 ) =PP*RR**3/C0N+0 . 5 *UU* *2 

EE (2)=(-PP*RR* *2*RPR+RR* * 3 * PPR ) /CON+UU*UPR 

EE (3) =(PP*RR*RPR**2-RR**2*PPR*RPRl /CON+O . 5*UPR**2 

EE(4)=(-PP*RPR**3+RR*PPR*RPR**2 ) /CON 

EE ( 5 ) =-PPR*RPR* * 3 /CON 
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c 
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RETURN 

END 


SUBROUTINE GAUSS (NINT , SAMP , WEIGHT) 

C 

COMMON/PARAMT/SOUND , GAMMA , ADM I , ADMO 
COMMON/PVALUE/XLENG , TLENG , IPERT , PI , EPSI , VISCO 
C 

DIMENSION SAMP(NINT) , WEIGHT(NINT) 

C 

N =N I NT 
M=(N+l)/2 
E1=N* (N+l ) 

DO 1 I =1 , M 
T=(4*I-1) *PI/(4*N+2) 

XO= ( 1 .0— (I. 0—1 .0/N) / (8 .0*N*N) ) *COS(T) 

PKM1=1 . 0 
PK=XO 

DO 3 K=2 , N 
T1=X0*PK 

PKP1=T1-PKM1- ( T1-PKM1 ) /K.+T1 
PKM1=PK 
3 PK=PKP1 

DEN=1 . -XO*XO 
D1=N* ( PKM1-X0* PX) 

DPN=D1/DEN 

D2PN= ( 2 . 0*X0*DPN-E1*PK.) /DEN 

D3PN= ( 4 . 0*X0*D2PN+ ( 2 . 0-El ) *DPN) /DEN 

D4PN=(6 . 0*X0*D3PN+ ( 6 . 0-El) *D2PN) /DEN 

U=PK/DPN 

V=D2PN/DPN 

CH=-U*(1.0+0.5*U*(V+U*(V*V-D3PN/(3.0*DPN) ) ) ) 

P=PK+CH* ( DPN+O . 5 *CH* (D2PN+CH/3 . 0* (D3PN+0 . 25*CH*D4PN) ) ) 
DP=DPN+CH* ( D2PN+0 . 5 *CH* ( D3PN+CH*D4PN/ 3 . 0) ) 

CH=CH-P/DP 
SAMP ( I ) =XO+CH 

CFX=D1-CH*E1* ( PK+O . 5*CH* (DPN+CH/3 . 0* (D2PN+0 . 25*CH* 

& (D3PN+0. 2*CH*D4PN) ) ) ) 

1 WEIGHT ( I ) =2 . * ( 1 . 0-SAMP ( I ) *SAMP ( I ) )/(CFX*CFX) 

MM=N/2 

DO 25 J =1 , MM 
IF(2*M.EQ.N) GOTO 22 
SAMP ( M+J ) =-SAMP ( M- J ) 

WEIGHT (M+J ) =WEIGHT ( M-J ) 

GOTO 25 

22 SAMP (M+J )=-SAMP(M+J.-J) 

^ WEIGHT (M+J ) =WEIGHT (M+l-J ) 

25 CONTINUE 

IF (M+M . GT . N ) SAMP (M) =0 . 0 

C 

RETURN 

END 

C 

c ; •• 

SUBROUTINE SHAPE (NINT, XI , PHI) 

C 

COMMON / PARAMT/ SOUND , GAMMA , ADM I , ADMO 


COMMON/PVALUE/XLENG , TLENG , I PERT , PI , EPSI , VISCO 
C 

DIMENSION XI(NINT) ,PHI(2,NINT) 

C 

DO 100 1=1 , NINT 

PHI (1 , I ) =0 . 5* (1 . 0-XI ( I) ) 

PHI (2, I)=0. 5* (1. 0+XI ( I) ) 

100 CONTINUE 
C 

RETURN 

END 

C 

C 

SUBROUTINE PRIME ( PZERO . UZERO , XX , TT , PPRIME , UPRIME) 

C 

COMMON / P ARAMT /SOUND , GAMMA , ADM I , ADMO 
COMMON/PVALUE/XLENG , TLENG , I PERT . PI , EPSI , VISCO 

C 

PPRIME=0 . 0 
UPRIME=0 . 0 

C 

DO 100 1=1 , 12 
CON=I*PI 

CN=2 . 0*SQRT ( 1 . O+CON : *2) /CON* *2 
PHASE=-ATAN ( 1 . O/CON > 

CKN=CON/XLENG 
OMEGA=CON * SOUND / XLEMG 
XCON=CKN*XX 
TCON=OMEGA*TT-PHASE 

PPRIME=PPRIME+CN*COo ( XCON ) * SIN ( TCON ) 
UPRIME=UPRIME+CN*SIis (XCON) *COS (TCON) 

100 CONTINUE 
C 

PPRIME=EPSI*PZERO*PFRIME 

UPRIME=EPSI*SOUND*UPRIME 

C 

RETURN 

END 


